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ABSTRACT 


An  efficient  finite  element  model  and  solution  technique  have  been 
developed  for  the  analysis  of  unrestrained  flexible  structures  undergoing 
large  elastic  deformations  coupled  with  gross  nonsteady  translational  and 
rotational  motions  with  respect  to  an  inertial  reference  frame.  The  nonli¬ 
near  coupled  differential  equations  resulting  from  the  finite  element 
approximation  are  integrated  timewise  using  an  implicit-explicit  split 
operator  numerical  integration  scheme  which  treats  the  stability  sensitive 
terms  of  the  equation  implicitly  while  the  rest  of  the  equation  is  treated 
explicitly.  The  motion  of  simple  spacecraft  structures  consisting  of 
flexible  beams  attached  to  rigid  masses  and  including  the  effect  of  control 
forces  has  been  studied  using  three-node  eighteen-degree-of-freedom  three 
dimensional  beam  elements  based  on  the  total  Lagrangian  description. 
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The  velocity  vector  is  defined  as  follows: 


11m  4 

At+0  At 


(2.11) 


Substituting  eq.  (2.10)  into  eq.  (2.11), 

Aq 


u_  =  1  im 

-  a 

(2.12) 

At  -►0 

But 

• 

q  =  lim 
At+0 

ia, 

At 

(2.13) 

Thus 

•  • 

u  =  N*  q 

(2.14) 

Similarly, 

V  =  N*  q+  N*  q  *  N*  q_  (2.15) 

•  • 

noting  that  N*  q_is  either  zero  for  solid  elements  or  very  small  in  the 
case  of  beam,  plate,  or  shell  elements. 


2.4  Conservation  of  Linear  Momentum 


If  Fj  is  the  sum  total  of  applied  forces  on  the  body,  then  conser¬ 
vation  of  linear  momentum  requires  that  the  following  equation  be 
satisfied: 

dR 


a  f 

at  Jp 


d2R 


p  Clt  av  =  1  p — “  av  =  -T  (2.16) 

v  v 

In  eq.  (2.16),  the  integral  is  defined  over  the  original  undeformed  con¬ 
figuration.  Substituting  the  expression  for  acceleration  eq.  (2.8)  into  eq. 
(2.16)  and  using  eqs.  (2.14)  and  (2.15),  eq.  (2.16)  becomes 


M(V0  ♦  n  V  )  +fTi+G_£  +  2  £  £T  i  +  Ic  =f.y 


(2.17) 


•.  .  -l 


A 


7 


0 


where  A  = 


0  -A, 


(2.7b) 


Applying  equations  (2.7)  to  equation  (2.6),  the  cross  products  can  be 
replaced  as  follows: 


dzR  *  —  •• 

— “  =V/>  +  nV_  +  u  +  nr  +  2fiu+nr 
dt 2  -o  --o  -  —  ~  - 


(2.8) 


2.3  finite  Element  Expression  for  Velocity  and  Acceleration 

In  the  finite  element  method,  the  structure  being  analyzed  is  divided 
into  a  finite  number  of  sections  or  elements.  Within  each  element  the 
exact  displacement  u_i$  approximated  by  polynomials  containing  unknown 
constants  which  generally  represent  the  displacement  at  a  finite  number  of 
points  or  nodes  within  the  element  and  on  its  boundaries.  It  is  these 
nodal  displacements  which  are  solved  for.  For  time  dependent  problems, 
the  displacement  u  can  be  written  in  an  incremental  form  as 


u  =  U  +  AU 


(2.9) 


In  equation  (2.9),  ‘'u  represents  the  value  of  _u_  at  time  t  while  A£  is  the 
change  in  displacement  between  time  t  and  t  +  At.  In  the  finite  element 
formulation  the  value  of  au  within  each  element  can  be  expressed  as 


Au_  =  N_*  Acj_ 


(2.10) 


Here  Aq  is  the  nodal  displacement  vector  and  IN*  is  the  matrix  of  shape 
functions  and  is  a  function  of  a  set  of  local  coordinates  within  each  ele¬ 
ment  as  well  as  the  initial  angular  displacements.  The  determination  of  N* 
is  given  in  Chapter  III. 


6 


In  order  to  express  the  motion  in  terms  of  the  rotating  body  axes,  it 
is  necessary  to  express  the  acceleration  of  point  P  in  terms  of  body 
axis  coordinates.  Thus  the  following  expression  for  the  acceleration  is 
used: 


Wo 

dt 


+r+nxr+2nxr+nx 


(o  *  I) 


(2.2) 


dV 

where  =  the  acceleration  of  the  origin  of  the  body  axis 

r  =  the  velocity  of  the  material  point  p  relative  to  the  body  axis 

V  =  the  acceleration  of  the  material  point  P  relative  to  the 

body  axis 

C  x  r  =  the  "tangential"  acceleration 
2n.x  £  =  the  Coriolis  acceleration 
£x  (n  x  r)  =  the  centripetal  acceleration 
In  terms  of  the  body  axis  coordinates. 


fi  x  V 
—  — o 


(2.3) 


Also  note  that  since  rQ  =  rQ  =  0, 

r  =  u.  (2.4) 

••  ••  .  . 

£  =  u_  (2.5) 

Substituting  equations  (2.3)  -  (2.5)  into  equation  (2.2), 

a2u  *  ...  • 

— ?  sV  +£xVn  +  u  +  fixr  +  2fixu  +  nx(fixr)  (2.6) 

dt2  "°  ~  ~  “  “ 

In  order  to  simplify  equation  (2.6)  a  bit,  the  cross  products  can  be 

replaced  by  matrix  products.  If  A  and  are  any  3-vectors,  then 

A  x  B  =  A  B  (2.7) 
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CHAPTER  II 


EQUATIONS  OF  MOTION 


2.1  Introduction 

For  elastic  structures  undergoing  gross  translational  and  angular 
motion  as  well  as  small  or  large  elastic  deformations  the  motion  can  be 
described  by  the  following  three  sets  of  equations: 

1.  The  conservation  of  linear  momentum  which  is  a  vector  equation 
describing  the  gross  translational  motion. 

2.  The  conservation  of  angular  momentum  which  is  a  vector  equation 
describing  the  gross  rotational  motion. 

3.  The  principle  of  virtual  work  which  describes  the  deformations. 

If  the  structure  is  very  flexible,  the  deformed  configuration  may  be  quite 
different  from  the  original  undeformed  configuration.  Thus  the  elastic 
deformations  will  be  coupled  with  the  gross  translations  and  rotations, 
especially  if  the  applied  loads  are  deformation  or  velocity  dependent. 

2.2  Geometry  and  Kinematics 

Consider  the  deformable  body  pictured  in  Fig.  1.  A  set  of  mutually 
orthogonal  axes,  x^ ,  x^,  and  x3  are  fixed  in  the  undeformed  body  at  point  0. 
Point  0  is  located  a  distance  RQ  from  a  set  of  mutually  orthogonal  inertial 
axes  Xj ,  X2,  X3,  centered  at  point  C.  The  body  axes  are  translating  with 
velocity  VQ  and  rotating  with  angular  velocity  £  relative  to  the  inertial 
axes.  A  point  P  located  a  distance  rQ  from  0  displaces  by  £  to  point  P'  as 
the  body  deforms  so  that  it  is  now  a  distance  £  from  0.  Point  P’  is 
located  a  distance  £  from  C  where  R  is  as  follows: 

i=i*o+I=!io+-o  +  -  (2,1) 
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flexible  structures  which  are  undergoing  large  elastic  deformations  coupled 
with  gross  nonsteady  rigid  body  translational  and  rotational  motions  with 
respect  to  an  inertial  reference. 

The  formation  and  solution  scheme  used  for  this  research  can  be 
briefly  described  as  follows.  The  governing  equations  of  motion  are 
derived  using  momentum  conservation  principles  and  the  principle  of  virtual 
work.  The  finite  element  approximation  is  applied  to  the  equations  of 
motion  and  a  matrix  form  of  those  equations  is  obtained.  The  resulting  set 
of  second  order  matrix  differential  equations  is  solved  timewise  by  direct 
numerical  time  integration  using  an  implicit-explicit  split  operator 
scheme.  This  scheme  treats  the  terms  which  control  the  stability  of  the 
solution  implicitly  while  the  terms  which  are  less  sensitive  to  stability 
are  treated  explicitly.  The  solution  technique  developed  is  tested  on 
simple  spacecraft  consisting  of  long  slender  uniform  beams  attached  to  a 
rigid  mass  and  modeled  by  three  dimensional  beam  elements.  The  effects  of 
control  forces  on  the  motion  of  the  spacecraft  are  also  considered. 


not  the  associated  modal  generalized  forces  vanish  identically.  McDonough 
[3]  considered  the  formulation  of  the  global  equations  of  motion  of  an 
unrestrained  deformable  body  using  translating  and  rotating  reference 
frames.  The  motion  of  the  body  (including  deformations)  is  of  unrestricted 
magnitude  in  the  analysis.  Fraejis  de  Veubeke  [4j  considered  the  motion  of 
a  flexible  body  undergoing  arbitrarily  large  rotations  with  respect  to  an 
inertial  frame.  The  motion  was  split  into  a  mean  rigid  body  motion  and  a 
relative  motion  taking  into  account  the  deformtions.  This  mean  rigid  body 
motion  is  chosen  so  as  to  minimize  the  mean  square  of  relative  displace¬ 
ments.  Kane  and  Levinson  [5,6)  considered  different  methods  for  for¬ 
mulating  the  equations  of  motion  for  complex  flexible  spacecraft.  These 
methods  included  momentum  principles,  D'Alembert's  Principle,  and 
Lagrange's  equations,  among  others.  They  also  developed  an  algorithm  for 
producing  numerical  simulations  of  large  motions  of  a  nonuniform  flexible 
cantilever  beam  in  orbit  using  the  finite  element  method.  Santini  [7J  has 
studied  the  stability  of  both  nonspinning  and  spinning  flexible  spacecraft 
in  a  gravitational  field  by  the  superposition  of  a  rigid  motion  plus  a  com¬ 
bination  of  structural  modes.  Other  investigtors  [8-17 J  have  also  studied 
these  types  of  problems  using  both  general  analyses  or  in  connection  with 
more  specific  types  of  flexible  spacecraft.  It  appears  however  that,  in 
spite  of  the  progress  made  in  the  analysis  of  such  problems,  little 
attention  has  been  paid  to  the  development  of  the  finite  element  method 
for  the  dynamics  of  unrestrained  structures  undergoing  large  elastic  defor¬ 
mations  coupled  with  nonsteady  gross  translational  and  rotational  motions. 

With  these  problems  in  mind,  the  objective  of  this  research  is  to  use 
the  finite  element  method  to  determine  the  time  response  of  unrestrained 


CHAPTER  I 


INTRODUCTION 


In  order  to  predict  the  motion  of  many  types  of  flexible  spacecraft, 
it  is  necessary  to  accurately  simulate  the  time  response  of  an  unrestrained 
structure  which  is  undergoing  large  elastic  deformations  as  well  as  gross 
nonsteady  rigid  body  translations  and  rotations.  For  such  structures,  the 
large  elastic  deformations  are  coupled  with  the  rigid  body  motions 
resulting  in  a  complicated  set  of  nonlinear  differential  equations.  Such 
spacecraft  may  be  simple  enough  to  be  modeled  as  rigid  bodies  supporting 
flexible  beams  or  they  may  be  more  complicated  structures  consisting  of 
frames,  plates  and  shells  in  combination  with  one  or  more  rigid  bodies. 

For  example,  consider  a  large  space  structure  consisting  of  a  frame  made  up 
of  long,  slender,  flexible  beams  connected  to  one  or  more  rigid  masses.  If 
such  a  spacecraft  were  to  execute  a  sudden  rotational  maneuver  or  reorien¬ 
tation,  then  large  elastic  deformations  coupled  with  rigid  body  motion 
would  occur.  An  accurate  time  response  analysis  of  the  motion  of  the 
structure  would  be  necessary  in  order  to  predict  the  orientation  of  the 
structure  especially  if  it  were  necessary  to  determine  the  pointing 
accuracy  of  any  sensors  which  may  be  attached  to  the  spacecraft. 

Extensive  research  has  been  done  in  the  field  of  the  dynamics  of 
flexible  spacecraft.  The  motion  of  unrestrained  flexible  structures  has 
been  discussed  by  Bisplinghoff  and  Ashley  [l]  who  considered  small  vibra¬ 
tions  of  aircraft  structures  using  a  modal  technique.  Ashley  [2]  also 
studied  gravitational  excitation  of  very  simple  elastic  spacecraft  under 
the  restriction  of  infinitesimal  elastic  displacements  as  well  as  cate¬ 
gorizing  typical  free-free  structural  configurations  according  to  whether  or 
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where 


p  dv  =  Total  mass  of  body 


(2.18) 


T  N 
PT  =  i 


p  N*  dv 


(2.19) 


6  =  -  p  r  dv 


(2.20) 


fc  =  n2  /  p  £  dv 


(2.21) 


The  nodal  accelerations  q  are  coupled  with  the  translation  of  the  body 
axes  and  the  £  matrix  represents  the  mass  of  the  finite  elements  in  the 
coupling  terms.  Note  that _P  is  the  summation  or  assembly  over  all  N  finite 
elements.  The  elements  of  the  £  matrix  represent  the  location  of  the 
instantaneous  center  of  mass.  The  vector  is  the  force  due  to  the 

centripetal  acceleration.  Details  of  the  P,  G,  and  f  matrices  are  given 
in  Appendix  C.  Rewriting  eq.  (2.18)  with  all  the  acceleration  terms 
grouped  together  results  in  the  following: 


where 


T  ••  T  #  * 

p  a  +  £  £  +  m  v0  =  fT 


It  1  It  -  M  '  2  “  a  -  lc 


(2.22) 


(2.23) 


2.5  Conservation  of  Angular  Momentum 


If  My  is  the  sum  total  of  the  applied  moments  about  the  inertial  axes, 
then  the  conservation  of  angular  momentum  requires  that  the  following 
equation  be  satisfied: 


*f>  "* 


p(R  x  dv  =  Mt 
dt2  1 


(2.24) 


Substituting  eq.  (2.1)  into  eq.  (2.24)  and  using  eq.  (2.16),  eq.  (2.24) 
becomes 


%  -  B0  x  Fj  +  J  p(r  x  ~)  dv 


(2.25) 


But  the  applied  moment  about  the  inertal  axes  can  also  be  written  in  the 
following  form: 


-  %  +  R0  x  It 


(2.26) 


Here  MQ  represents  the  applied  moment  about  the  body  axes.  Thus  equation 
(2.25)  reduces  to 


M0  -  /  p  (r  x  — )  dv 


(2.27) 


Now,  substituting  eq.  (2.8)  for  the  accelerations  into  eq.  (2.27),  and 
using  eqs.  (2.18)  -  (2.21),  eq.  (2.27)  becomes  as  follows: 


Uo  ■  §Ti0  *  aT  zio  *  aT  a*  It  a*  aCor  *  a 


(2.28) 


where 


IT  3  -  /  P  r2  dv 


(2.29) 


N  f  _ 

i  /  p  r  N*  dv 

1-1  J 


(2.30) 
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(2.31) 


M.rtr  =  2 
—cor 


f 

/  p  £  n_  u_ 


dv 


"cent  -  /  p  r  n2  x  dv 


(2.32) 


The  matrix  represents  the  instantaneous  moment  of  inertia  matrix  of 

the  entire  structure.  The  hJ  matrix  represents  the  distributed  moment  of 

inertia  matrix  of  the  finite  elements  and  when  multipled  by  the  nodal  acce- 
•  • 

leration  vector  q  provides  the  coupling  term  between  the  rotary  inertia  due 
to  the  elastic  displacements  and  the  rotation  of  the  body  axes.  The  M... 

"“vUI 

and  *  vectors  represent  moments  due  to  the  Coriolis  and  centrifugal 
forces  respectively.  Details  of  £y,  H,  M(.or  and  Mcent  are  provided  in 
Appendix  C.  Expressing  eq.  (2.28)  in  a  form  similar  to  eq.  (2.23)  results 
in  the  following: 


ht  i  +  It  £  +  £  \  = 


where 


m,  =  M  -  G  V  -  M  -  M  . 
—i  -o - o  —cor  -cent 


(2.33) 

(2.34) 


2.6  Principle  of  Virtual  Work 

The  principle  of  virtual  work  is  used  to  describe  the  elastic  displa¬ 
cements  of  the  structure.  Using  tensor  notation,  the  principle  of  virtual 
work  for  a  solid  body  undergoing  large  deformations  can  be  written  as 
follows: 


6u  ds  *  0 


(2.35) 


10 


„• 


The  first  Integral  represents  the  virtual  work  done  by  the  internal  forces 
where  S  is  the  second  Pyola-Kirchhoff  stress  tensor  and  6E_  is  the  virtual 
Green  strain  tensor.  The  second  integral  represents  the  virtual  work  done 
by  the  body  force  The  virtual  displacement  is  represented  by  6u.  The 

third  integral  represents  the  virtual  work  done  by  the  applied  surface 
tractions  T.  The  portion  of  the  surface  over  which  the  tractions  are 
applied  is  S  . 

Now,  consider  the  second  integral  in  eq.  (2.35).  The  body  forces 
acting  on  the  body  are  as  follows: 


(2.36) 


The  first  term  represents  the  inertia  force  while  £b  represents  any  other 
applied  body  forces  such  as  gravity.  If  eq.  (2.36)  is  substituted  into  eq. 
(2.35)  and  if  eqs.  (2.8),  (2.10),  (2.14),  (2.15),  (2.18)  -  (2.21),  and 
(2.29)  -  (2.32)  are  used,  then  the  second  integral  in  eq.  (2.35)  becomes  as 
follows  for  a  single  element: 


F  •  6udv  =  «£T  (M  i  +  H  £  +  _P  Vq  +  C  £  +  P  £  VQ  +  Fc  -  fB)  (2.37) 


where 


M  =  J  p  N* 1  N*  dv 
v 

C  =  2  f  p  N*T  n  N*  dv 


(2.38) 

(2.39) 


Ij.  =  I  p  £2  (£0  +  u)  dv 


(2.40) 


» 


t 
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(2.41) 


-B 


Jf T  Fb  dv 


Note  that  _M  is  the  usual  consistent  mass  matrix.  The  £  matrix  is  the 
gyroscopic  matrix  which  is  skew- symmetric  and  represents  the  contribution 
of  the  Coriolis  acceleration  to  the  inertia  force.  The  Fr  vector  is  the 
force  vector  due  to  the  centrifugal  force.  Finally,  Fg  is  the  applied  force 
term  due  to  gravity.  Details  of  eqs.  (2.38)-(2.41 )  are  given  in  Appendix 
D.  In  addition,  the  volume  integrals  are  evaluated  over  the  original  unde¬ 
formed  volume. 

The  third  integral  in  eq.  (2.35)  is  treated  in  the  standard  fashion 
for  finite  elements. 


T 


ds  = 


lo 


(2.42) 


Here  F^  is  the  vector  of  applied  forces  at  the  nodes  which  result  from  the 
applied  tractions. 


2.6.1  The  Elastic  Stiffness  Matrix 

The  application  of  the  finite  element  approximation  and  a  trapezoidal 
rule  integration  to  the  first  integral  in  the  principle  of  virtual  work 
leads  to  a  system  of  nonlinear  equations  that  can  be  solved  by 
Newton-Raphson  iteration.  With  this  iteration  in  mind,  the  displacements 
the  second  Pyola-Kirchhoff  stress  and  the  Green  strain  tensor  £  are 
written  in  incremental  form  as 


o 

u  =  _u  + 

AU 

(2.43a) 

+ 

cT' 

11 

AS 

(2.43b) 

E  =  °£  + 

AI 

(2.43c) 
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where  u  =  Displacement  vector  at  state  i 

°S  =  Second  Pyola-Kirchoff  stress  vector  at  state  i 
E  =  Green  strain  vector  at  state  i 

Au  =  Incremental  change  in  u_  between  states  i  and  i  +  1 

AS  =  Incremental  change  in  £  between  states  i  and  i  +  1 

aE  =  Incremental  change  in  £  between  states  i  and  i  +  1 

Note  that  i  refers  to  the  ith  iteration  in  the  Newton-Raphson  iteration. 

Substituting  eq.  (2.43)  into  the  first  integral  in  eq.  (2.35), 


5ET  S  dV  =  /  6  ET  °S  dV  +  /  6  ET  AS  d V 


The  constitutive  law  may  be  written  as 


(2.44) 


AS  =  C  AE 


(2.45) 


where  C  =  Material  property  matrix 


The  material  property  matrix  C  may  in  general  depend  upon  state  n.  All  the 
materials  considered  here,  however,  are  linearly  elastic  which  results  in 
the  £  matrix  being  constant. 

Substituting  eq.  (2.45)  into  eq.  (2.44), 


6  ET  S  dV  =  /  6  ET°S  dV  +  /  <5  ET  C  AE  dV 


(2.46) 


Now,  using  tensor  notation,  the  Green  strain  tensor  can  be  written  as: 


,  3u.  3u .  au.  au. 

F  -  ( L  +  J  +  *  *  \ 

tij  ~  7  '"5X7  axi  axi  aXj' 


|  -j 

also  6Ei4  =  i  (-w— 


asu.  afiu .  au.  asu.  au.  asu. 

_ 1  X  _ i  X  _ 1  i  X  _ 1  _ 1' 


i  j  =  7  +  axi  +  axi  axi  +  ax .  ax . 

J  •  *  J  J  ^ 


(2.47) 


(2.48) 


v>:> 
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and  u^  =  u^  +  au^ 


(2.49 


Substituting  equation  (2.49)  into  equations  (2.47)  and  (2.48)  results  in 


Eij  *  E>j  *  SEu 


8Eij  =  6tij  +  '"ij 


.  3°u,  sV  jV  >\ 

where  E^.  =  7  (-^-  +  +  -55 r  -gH 

J  *  I  J 


Aeij  +  Anij 


,  3AU.  3AU.  3  U.  3AU.  3  U.  3AU. 

i _ 1  + _ J  + _ £ _ £  + _ £ _ £\ 

2  v  aXj  3Xi  3Xi  3Xj  3Xj  3Xi  ’ 


.  36u.  35U.  3°u.  35u,  3°u,  35u. 

.  .  =  J.  f  '  +  _ i  +  - _  —  _  -  _  +  - -  - .  \ 

1J  2  3Xi  3Xi  3Xj  3Xj  SX^ 


T  3AU  3AU 

.  T  k  k 


i  j  “  2  ”57"  ~57~ 


.  3AU  36U.  3AU.  35U. 

6"ij 


(2.50) 


(2.51) 


(2.52) 


(2.53) 


(2.54) 


(2.55) 


(2.56) 


(2.57) 


Note  that  An^-  is  quadratic  in  au^.  Now,  substituting  equations  (2.51)  and 
(2.53)  into  equation  (2.46)  and  using  matrix  notation  results  in 


6E  S  dv 


-  ( 


6e"^  C  Ac  dv  +  f  5rJ°S  dv  +  f  6e^°S  dV  +  A„ 


where  An  =  Higher  order  terms,  i.e.,  terms  nonlinear  in  aik 

2.6.2  Finite  Element  Approximation 

The  finite  element  approximation  to  the  incremental  displacements  bet 
ween  iteration  steps  i  and  i  +  1  is  as  follows: 


Au  =  Aq 

where  =  Shape  function  matrix  dependent  upon  state  i 

Aq  =  Nodal  displacement  between  states  i  and  i  +  1 


(2.59 


The  virtual  displacement  6u  can  thus  be  expressed  as  follows: 


6u  =  N*  fiq 


(2.60 


Also,  as  will  be  shown  in  Chapter  3,  it  is  possible  to  express  the  incre¬ 
mental  strains  Ae  within  each  element  as  follows: 


Ae  *  B 


(2.61 


Also,  fie  =  B  fiq 


(2.62 


where  Ae^  =  LAe^  AE22  Ae^j  2Ae^  2Ae^ 


(2.63 


In  addition,  it  will  be  shown  that 


finToS  =  fiaJ  Bs  A£ 


(2.64 


If  equations  (2.61)  -  (2.64)  are  substituted  into  the  integrals  on  the 
right-hand  side  of  equation  (2.58),  these  integrals,  when  taken  over  the 
volume  of  the  ith  element  become: 


f  fieT  C  Ae  dv  =  fig1  (  f  B T  C  B  dv)  A£  =  6^  Ki  Ag, 


(2.65 


J  «n  S 


dv  =  fig'  (I  Bs  dv)  Aq  =  6q  K  acj. 


J  6  eToS  dv  =  fig1  ( J  BT°S  dv)  =  fig_T 


where  K-  =  /  BT  C  B  dV  =  element  basic  stiffness  matrix 


K.  =  /  B„  dV  =  element  initial  stress  stiffness  matrix 

-si  /  ~s 


F„  =  /  BToS  dV  =  element  initial  stress  force  vector 

-Si  -  - 


Thus,  if  the  higher  order  terms  Ap  are  neglected,  the  virtual  work 
expression,  equation  (2,58),  becomes  when  summed  over  all  N  elements 


6ETS  dv 


i=E1  6ai[(!ii  +  !i5i)  a9i  +  fs.] 


Assembling  the  elements  in  the  usual  way. 


f  5E^S  dv  =  fi£^  (K^  A9  +  F^) 


where 


(2.73) 


s  -  +  -s  =  Total  9^obal  elastic  stiffness  matrix 


F<.  =  Global  initial  stress  force  vector 


(2.74) 


Finally,  by  using  eqs.  (2.71),  (2.42),  and  (2.37),  the  virtual  work 
expression  can  now  be  written  as  follows: 


+  Llo  +  L  i  +  JSe  A±=Lq 


(2.75) 


where  f q  "  lo  +  £b  -  fc  -  Is  -  f  SL  1q 


(2.76) 


Thus,  Eq.  (2.75),  (2.33)  and  (2.22)  are  the  three  sets  of  nonlinear  coupled 
equations  of  motion  for  the  unrestrained  system  shown  in  Fig.  1. 


The  next  section  developes  the  element  stiffness  matrices  and  initial 
stress  force  vector  for  the  three  node,  18  degree  of  freedom  beam  element. 
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CHAPTER  III 


DEVELOPMENT  OF  THE  STIFFNESS  MATRIX  FOR  THE  THREE  NODE  BEAM  ELEMENT  “• 

3.1  The  Three  Node  18  Degree-of-Freedom  Beam  Element  Geometry 

In  Figure  2  is  shown  a  section  of  the  beam  both  before  and  after 
deformation.  It  is  assumed  that  during  the  large  deflection  that  there  can 
be  large  rotations  but  that  the  strains  remain  small.  Also,  it  is  assumed 
that  there  is  no  warping  of  the  cross  section,  i.e.,  plane  sections  remain 
plane.  The  following  quantities  are  then  defined: 


»  *2*  —3  "  a  set  mutua11y  orthogonal  axes  in  the  undeformed 
beam,  centered  at  point  0  with  a^  tangent  to  the 
reference  line  (beam  axis) 

aj,  a2’  -3  “  a  set  of  mutual  orthogonal  axes  in  the  deformed  beam 
centered  at  O'  with  a^  and  a^  lying  in  the  deformed 
cross-sectional  plane 
XQ  -  Position  vector  of  point  0 
X  =  XQ  +  rp  =  Position  vector  of  a  general  point  P  in 
the  plane  of  the  cross-section 
uQ  =  Displacement  of  point  0 
u  =  Displacement  of  point  P 

Note  that  both  a£  and  a^  as  well  as  a^  and  a^  remain  in  the  cross-sectional 
plane  because  of  the  assumption  that  plane  sections  remain  plane.  Thus  £] , 
a, and  represents  the  orientation  of  the  a_j ,  a^,  axes  after  their 
rotation  due  to  the  deformations. 

From  the  geometry  in  Figure  2, 
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Next,  define  a  local  coordinate  system,  x,  y,  z,  such  that  x,  y,  and  z  are 
the  a^ ,  a^,  and  a^  coordinates  respectively.  Thus, 

=  y  a2  +  z  a3  (3.2) 

and  A  *  £o  +  £p  (3.3) 

Also,  because  the  cross  section  only  translates  and  rotates  (and  doesn't 
distort) , 


=  y  +  2  ±2  (3.4) 

Substituting  equations  (3.2)  and  (3.4)  into  equation  (3.1)  yields, 

1  =  Ho  +  ‘  -2*  y  +  ^-3  '  ^3^  2  (3.5) 

In  the  Total  Lagrangian  description  all  variables  are  referred  to  the  ori¬ 
ginal  undeformed  configuration.  Thus,  it  is  necessary  to  determine  a£  and 
a3  in  terms  of  ^ ,  a2  and  a3.  In  order  to  do  that,  a  way  of  describing  the 
rotations  must  be  found.  There  are  many  well  known  ways  of  doing  this, 
probably  the  best  known  way  being  Euler  angles  (Fig.  3a).  Another  way  of 
doing  this  is  by  using  space-three  or  body-three  angles  (Figs.  3b,  3c) 

[34].  In  the  space-three  1-2-3  description,  a_j,  £2  and  _a3  are  rotated  suc¬ 
cessively  about  three  separate  axes  as  follows: 

1.  First,  rotate  by  an  amount  about  the  ori gi nal  ^  axis. 

2.  Next,  rotate  the  new  configuration  0£  about  the  original  a^  axis. 

3.  Finally,  rotate  the  last  configuration  @3  about  the  original  a3 
axi  s. 

The  body-three  1-2-3  sequence  is  similar  except  that  instead  of  rotating 
about  the  original  axes,  the  rotation  is  about  the  body  axes.  It  is  also 
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possible  to  use  sequences  other  than  1-2-3,  for  example,  a  space-three 
2-3-1,  or  body-three  1-3-2  or  any  sequence  which  results  in  the  required 
orientation  is  possible.  In  the  case  of  a  1-2-1  or  3-1-3  sequence  or  any 
sequence  in  which  there  is  a  rotation  about  only  two  axes,  the  rotation 
would  be  referred  to  as  a  space-two  or  body-two  rotation.  Note  that  the 
Euler  angles  correspond  to  a  body-two  3-1-3  rotation. 

For  the  purpose  of  this  research,  the  space-three  1-2-3  rotation 
sequence  was  chosen  to  describe  the  rotations  because  it  seemed  the  most 
straight  forward  and  because  none  of  the  others  offered  any  particular 
advantages  over  it. 

Once  a  method  for  describing  the  rotations  is  chosen,  the  relationship 
between  aj ,  a^,  a^  and  a^t  a 2,  a^  Cdn  be  determined.  For  the  space-three 
1-2-3,  this  is  as  follows: 


aj  =  cose2  cose3  a^  +  cose2  sine3  a2  -  sine2  a^  (3.6a) 

al,  =  (sinej  sine2  cose3  -  sine3  cose^)  £j  +  (sine1  sine2  sine3 


+  cose3  cose^)  £2  +  s i n cose2  a3 


(3.6b) 


£3  =  (cose^  sine2  cose3  +  sine3  sine^)  a^  +  (cose^  sine2  sine3 

-  cose3  sinOj)  a2  +  cos9|  cos02  a3  (3.6c) 

Next,  the  displacements  and  rotations  are  written  in  incremental  form 

u  =  °u  +  au_  (3.7) 

“o  =  \  +  ^0  (3*8) 

eR  =  °ek  +  a 9r,  k  =  1,2,3  (3.9) 
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-■f— *_■_  *_• . 


•  1 


■v 


T'  V 


If  equations  (3.6),  (3.8)  and  (3.9)  are  substituted  into  equation  (3.5) 
and  higher  order  terms  are  neglected  then  the  following  results: 


°u  =  0uQ  +  y  (°a^  -  a2)  +  z  (°a^  -  a3) 


=  &id0  +  ^  -  +  z 

where  °a!  =  a !  ("e^  °02,  °e3)  i  =  1,2,3 


(3.10) 

(3.11) 


Aj)  =  LA01  A02  A03J 
D  =  D  (°0lf  °02,  °03) 

E  =  E  (°0lt  °e2,  °03) 

The  elements  of  the  3x3  matrices  D  and  E  as  well  as  expressions  for 
°a!  -  ai ,  i  =  2,3  are  given  in  Appendix  A. 

3 .2  Finite  Element  Approximation 

The  three  node  eighteen  degree  of  freedom  element  is  pictured  in 
Figure  4.  The  local  coordinates  £,  y,  z  are  normalized  such  that  the  ori¬ 
gin  is  at  the  center  and  the  5  coordinate  varies  between  -1  at  one  end  of 
the  beam  and  +1  at  the  other  end.  The  transformation  between  the  global 
coordinates  and  the  local  coordinates  is  thus  as  follows: 


i*„E,  MOX^yfe^i, 


(3.12) 


Here  are  the  coordinates  of  the  three  nodes  and  N^(?)  are  the  shape 
functions  which  are 


1 


i^U)  *  J  c  u  -  u 

n2U)  =  1  -  c2 


(3.13a) 

(3.13b) 
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Thus  the  integral  is  first  evaluated  over  the  original  cross-sectional  area 


of  the  undeformed  beam  element  and  then  along  its  length  Next,  eq. 
(3.83)  is  substituted  into  eq.  (3.85)  and  the  various  matrix  products  and 
area  integrals  carried  out.  It  is  assumed  that  the  longitudinal  axis  of  the 
beam  is  through  its  centroid  and  that  the  cross  section  is  symmetric.  The 
following  equation  then  results: 


[A  (B£)T  C£  +  IZ2  vo,, ■  u 


(B f ) T  C£  B?  +  I 


yy  '—2 


(B/)T  c£  B£]  d i  (3.86) 


where 


(3.87) 


I 


zz 


y2  dA 


(3.88) 


(3.89) 


Note  that  A,  I22,  and  I  are  merely  the  area  and  area  moments  of  inertia 
of  the  beam  cross-section.  The  integration  along  the  length  of  the  beam  is 
done  numerically  using  Gaussian  integration  [28]: 


(3.90) 


The  a.  are  the  Gaussian  integration  points  and  H.  represents  the  weighting 

J  J 

factors  while  n  is  the  number  of  integration  points  used.  Since  the 
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For  an  isotropic  material,  the  local  3x3  matrix  of  material  properties  is 
then 


E 


C£ 


0 


0 


0  0 

BG  0 
0  BG 


(3.81) 


Here  E  is  the  Young's  modulus  of  the  material,  G  is  the  shear  modulus,  and 
B  is  the  shear  correction  factor.  For  an  isotropic  material  with  rec¬ 
tangular  cross  section,  B  is  equal  to  5/6. 


3.7  Calculation  of  the  Stiffness  Matrix  and  Initial  Stress  Force  Vector 
It  is  now  possible  to  calculate  the  element  stiffness  matrices  and 
initial  stress  force  vector  as  given  by  eqs.  (3.46)  -  (3.49).  First,  note 
that  the  local  B  matrices  are  formed  via  eqs.  (3.44)  as  follows: 

Bj  -  I  Bi*  i  =  0,1,2  (3.82) 


Thus  B_£  =  +  y  B*  +  z  B£  (3.83) 

Also  note  that  before  the _B$  matrices  can  be  formed,  it  is  necessary  to 
transform  the  stresses  from  the  local  coordinates  using  eq.  (3.40): 


°s  =  tt  V 


(3.84) 


The_PSi  matrices  can  then  be  formed  using  the  elements  of  the  vector 


from  which  in  turn  the  BSj  matrices  can  be  formed.  To  calculate  the  stiff¬ 


ness  matrix,  first  note  that  the  volume  integral  in  eq.  (3.46)  can  be  writ¬ 
ten  as  follows: 


(B£)T  C£  B£  dV  = 


(B£)T  C£  B£  dA  d£  (3.85) 
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along  its  length.  Thus,  in  a  manner  analogous  to  eq.  (3.70),  the  Sq 
vector  can  be  written  as 


also. 


°S  =  °s  +  y  °s  +  Z  °S 

do  J  i y  -Z 


fs=fs0  +  y^  +  z^ 


(3.77) 


(3.78) 


Finally,  if  eqs.  (3.59)  and  (3.78)  are  substituted  into  eq.  (3.76)  and  the 
matrix  products  carried  out,  the  following  expression  for  B  results: 


is  =  is0  +  y  Ssi  +  z  is2  +  yz  is3  +  y2  is 4  +  is 


(3.79) 


where  B  =  a'  Ps  An 

— Sq  —0  —  5q  —0 


Ssj  •  ij  £s0  io  +  il  -%  io  *  £  £s0 

*  il  !s0  io  *  il  £sz  A,  *  il  £s0  iz 

2s 3  ■  2j<Vz<£$2V  *  ^(£s0iz*£sA)  + 

is4  ■  4  £s0  iy  *  al  £sy  4,  +  £sy  io 

2s5  =  iz  £s0  iz  +  io  £sz  iz  +  iz  £sz-io 


(3.79a) 


(3.79b) 


(3.79c) 


(3.79d) 


(3.79e) 


(3.79f) 


3.6  The  Material  Property  Matrix 

The  material  property  matrix  C_  wi 1 1  in  general  depend  upon  the  state 
of  the  structure  and  as  such  it  is  a  function  of  time.  All  of  the 
materials  considered  here,  however,  are  linearly  elastic  which  results  in 
the  £  matrix  being  constant.  The  form  of  the  local  incremental  stress 
strain  relation  is  as  follows: 


AS*  =  C*  A e* 


(3.80) 
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(3.70) 


m,  u  *  ii  u  i  !■  i"  j'  *  * 


B  =  B0  +  y  B,  +  z  §2 

Next  it  is  necessary  to  determine  the  form  of  the  B$  matrix.  From  eq. 
(2.67), 

6nT°S  =  6£T  Bs  Aq_  (3.71) 


Using  the  expression  for  Sn  given  in  eq.  (2.57)  it  is  possible  to  show  that 


the  following  is 

true: 

6nToS  = 

,350,1 

1-57) 

i-s 

where 

£o 

0  0 

fs  = 

0 

p0  ° 

0 

£  Eo 

°s 

all 

°s  °s 

b12  b13 

Eo  = 

°s 

"*12 

°S  °s 

A22  A23 

°S 

n3 

°S  °S 

^23  A33 

Thus, 

using  eqs. 

(3.58) 

and  (3.63), 

6nToS  = 

%T  at 

Ps  A 

(3.72) 


(3.73) 


(3.74) 


or 


B$  =  A1  Ps  A 


(3.75) 

(3.76) 


Since  the  incremental  stresses  at  each  step  are  calculated  by  using  the 
constitutive  law  eq.  (2.45)  and  are  dependent  upon  the  strains  &e,  it 
follows  that  °S  will  vary  over  the  cross  section  of  the  beam  as  well  as 
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(3.66) 


V?)  (°a‘ 


S31  i 


Using  eqs.  (3.64)  -  (3.66)  as  well  as  eq.  (3.52),  the  elements  of  the 
expanded  De  matrix  can  be  obtained. 


D  =  D  +  y  D  +  z  D  (3.67) 

— e  — eo  ~€y  — ez 

The  terms  of  D^,  D  and  are  given  in  detail  in  Appendix  B. 

Now  we  are  prepared  to  calculate  the  matrix  given  in  eqs.  (2.61)  and 
(2.62).  Substitute  eqs.  (3.58),  (3.59),  and  (3.67)  into  eq.  (3.50)  and 
perform  the  matrix  multiplications  to  obtain  the  following: 


A_e  =  _B  acl 

where  j*  *  B0  +  y  -1  +  2  B2  +  •yz  -3  +  ^  -§4  +  -5 


So  ■  °e„  *„ 


B,  =  D  A  +  D  A 

-1  -Ey  -o  *e0  -y 

B?  =  D  A  +  Dr  A7 

—  c  — £2  — -O  — ‘Gq  """Z 


B.  =  D,  K  +  D  A, 
—o  -e2  ~y  -ey  —  Z 

*4  "  -ey  -y 

Be  =  D  A7 

—3  — €z  — Z 


(3.68) 

(3.69) 
(3.69a) 
(3.69b) 
(3.69c) 
(3.69d) 
(3.69e) 
(3.69f) 


Because  they  are  fairly  insignifient  and  add  unnecessarily  to  the 
complexity  of  the  stiffness  matrix,  the  higher  order  terms  yz  B3,  y2  and 
z2  Bg  are  neglected.  Thus  the  B  matrix  can  be  written  as  follows: 
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J-1 

ol 

-0 

j-1 

-0 

J-1 

2o 

(3.60) 


J"1 

D1 

— 

-y 

J’1 

Dv 

— 

-y 

J’1 

-y 

(3.61) 


j-1 

D1 

-z 

j-1 

£ 

1  '-3 1 

£ 

(3.62) 


It  should  also  be  noted  that  an  expression  similar  to  eq.  (3.58)  holds  for 

{iiJi}.  Thus, 

3X 


^  ■  A  6a 


(3.63) 


The  next  thing  needed  to  do  is  to  determine  the  elements  of  the  D£  matrix. 

.o 

This  will  require  the  evaluation  of  derivatives  like  — These  terms  can 
be  determined  by  first  differentiating  the  expression  for  °u  given  in  Eq. 
(3.14)  as  follows: 

ao  3  3N.  3  3N. 

Is  =  i=E,  -96  2«1  *  ,E,  —  fy  1  iz  -  iih  +  2  <  23  -  23>i^(3.64) 


N1«>  <°2z  *  — 2>i 


(3.65) 
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The  next  step  in  the  derivation  is  to  express  the  vector  of  derivatives 

terms  of  the  nodal  displacements  Aq.  This  can  be  done  by  differen¬ 
tiating  the  expression  for  the  finite  element  approximation  to  ^u  given  in 
eq.  (3.15). 

3AU 


3  3N-  3  3N. 

=  ij1  -5J  (yfii  +  zIi) 


35  i=l 


3AU 


3£ 


"S7*,f,  V<>£i 


3AU 


i=l 


az  if1  Ni^5^  -i  A-i 


By  using  eqs.  (3.53)  -  (3.55)  as  well  as  eq.  (3.24),  it  now  becomes 

3AU 

possible  to  write  as  follows: 


(3.53) 


(3.54) 


(3.55) 


3AU 


-1 


Aq 


(3.56) 


where  Dk  =  (5)  +  y  (5)  +  z  Dk  (5);  k  =  1,2,3  (3.57) 


The  elements  of  DK  are  composed  of  the  elements  from  the  and 

matrices,  as  well  as  terms  from  the  shape  functions  and  their  deri va¬ 
le 

tives.  The  exact  composition  of  the  D  matrix  is  given  in  appendix  B. 
Equation  (3.56)  can  be  rewritten  as  follows: 


3Au 

Hr*  =  £L A 1 


where  A  =  A  +  y  A ,  +  z  A, 
—  -0  — y  — z 


(3.58) 


(3.59) 
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(3.48 


£s<  * 


(B  )  S  dv 


3.49) 


3.5  Calculation  of  the  B  and  Bo  Matrices 

The  calculation  of  the  _B  and  matrices  used  in  the  integrals  of 
equations  (3.46)  -  (3.49)  is  now  considered.  The  incremental  strain  Aei  - 
given  in  eq.  (2.59),  when  written  in  vector  form  as  in  eq.  (2.66), 
can  be  expressed  as  the  following  matrix  product: 

3AU 

A-  =  -e  Hx*  <3 


3AU  3AUj  3AU1  3AU j  3AU2  3AU2  3AU2  3AU3  3AU.J  3AU3 

where  {_^=}  = 


3°Ui 

1  + 

3°U2 

3°U„ 

0 

0 

“31^ 

0 

0  ax' 

0 

3°u1 

3X2 

3°U? 

3°U~ 

0 

0 

0  1 

+  3X2 

0  0 

0 

3X2 

A, 

3°U? 

0 

0 

0 

0 

0  1 

3°u1 

3X2 

3°U1 

1*  1 
3X1 

0 

3°u2 
l+  n2 

3°u2 

3X} 

3°U. 

0 

3°u3 

3X1 

3X3 

0 

3°u. 

u-aq 

3°U2 

-5x7 

0 

_o  .0 

3  U,  3  U, 

c  1+  J 
3X1  3X3 

0 

0 

3°U1 

3X3 

3  °U-| 
3X2 

0 

3°U~  3°U? 

1  0 

3°U, 

1  + 

(3.52 
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where 


J 


The  elements  of  T  are  given  in  Appendix  A.  Since  all  the  other  quan¬ 
tities  involved  in  equation  (3.34)  are  tensors,  similar  expressions  hold 
for  them  also. 


Ae*  =  T  Ae 

(3.37) 

6e*  =  I  6e 

(3.38) 

fin*'  =  T_  6r^ 

(3.39) 

V  =  T°s 

(3.40) 

Next,  substitute  equations  (2.64)  - 

(2.66)  into  equations  (3.37)  -  (3.40) 

Ae4  =  B*  Aq 

(3.41) 

5e*  =  B*  6g_ 

(3.42) 

(6n£)TV  =  fiqT  85  Aq 

(3.43) 

where  B*  =  T  B 

(3.44) 

2s  ■  2S 

(3.45) 

Finally,  the  stiffness  matrices  and 

equations  (2.73)-(2.75) ,  become. 

the  initial  stress  force  vector. 

Ki  =  f  (B*)T  C*  B*  dv 

(3.46) 

3.4  Assumptions  on  Stress  and  the  Transformation  Matrix 


The  following  assumptions  have  been  made  concerning  the  local  stresses 
in  the  beam: 


322  I  I  S33  I  I  S23 


Here  the  i  superscript  refers  to  the  local  a^ ,  a^,  a^  coordinate  system. 

0  0  0 

The  S22.  S33,  and  S22  stresses  are  thus  considered  small  enough  to  be 
neglected. 

Since  the  virtual  work  integral  is  a  scalar,  it  is  invariant  under  any 
coordinate  transformation. 


6  E*T  S£  dv  =  /  6  ET  S  dv 


(3.34) 


The  transformation  from  local  to  global  coordinates  using  tensor  notation  is 
as  follows: 


Ers  "  (in  •  £r)(in  *  is)  Er 


(3.35) 


where  i  =  global  unit  vectors 


For  example. 


l  2  2  2 

E11  ‘  all  E11  +  a12  E22  +  a13  E33  +  2all  a12  E12 


+  2all  al 3  E1 3  +  2a12  a!3  E23 


Similar  expressions  hold  for  E^  and  e|3>  A  transformation  matrix  _T  can 
thus  be  defined  relating  the  local  and  global  strain  components: 


E£  =  T  E 


(3.36) 


V -.1 
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(3.28a) 


*■■  ■  p 


1  .  *r.  ~ 


ax 

H 


ax 

ay 


-2 


(3.28b) 


3X 

—  *  a.  (3.28c) 

az 

Thus,  the  Jacobian  matrix  is  as  follows: 


1  a 

7  311 

1  a 

7  a12 

£ 

7 

a13 

i  = 

a21 

a22 

a23 

*31 

a32 

a33 

(3.29) 


To  determine  J“^ ,  it  is  first  necessary  to  evaluate  the  determinate  of  the 
Jacobian.  First,  note  that 


1  i  1  =  7  £i  *  (^2  x  £3) 

(3.30) 

But 

£2  x  — 3  =  ~1 

(3.31) 

Thus 

1  i  1  =  1  £1  *  £1 

(3.32) 

and 

1  l\  =  7 

(3.33) 

This  expression  for  the  determinant  of  the  Jacobian  will  prove  useful  in  the 
evaulation  of  various  integrals. 
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where 


J=  Jocabian  matrix 


Solving  for  the  desired  derivative, 


f 

f  \ 

3 

3 

3X^ 

H 

3 

J 

3  \ 

iq 

3y  / 

3 

3 

32 

S.  J 

In  order  to  evaluate  the  Jacobian,  equation  (3.12)  is  used. 

*-  j,  "tUUot  ♦»%***3 


(3.24) 


(3.25) 


For  straight  beams,  however, 

*o2  -  X0,  ♦  {  il  (3.26a 

*03  "  *0|  *  '  it  (3.26b 

where  i  =  Length  of  the  beam  element 

Substituting  equations  (3.26)  into  (3.12)  and  using  the  expressions  for  the 
shape  functions  given  in  equations  (3.13)  results  in 

Is  lo,  +|  U+l)  a,  +  y  ^  +  z  a3  (3.27) 

Using  equation  (3.27)  to  evaluate  the  derivatives  in  the  Jacobian, 


AU  =  A£ 

where  N**NJ  +  yN*  +  zN* 

%  -  a  n2-3  a  n3i3] 

n*  =  Co  Ni21  o  n2d2  o  n3d3] 

n*  =  Co  N]E1  o  n2e2  o  n3e3] 


(3.17) 

(3.18) 

(3.19) 

(3.20) 

(3.21) 


Note  that  N*  N*  and  N*  are  3  x  18  matrices  since  U  is  the  3x3  iden- 

***»  "j  “"t  “O 

tity  matrix  and  0  is  the  3  x  3  null  matrix. 

3.3  Derivatives  and  the  Jacobian 

In  order  to  evaluate  the  integrals  in  the  virtual  work  expression,  it 
is  necessary  to  evaluate  terms  that  include  derivatives  with  respect  to  X_. 
Since  the  displacements  are  expressed  in  terms  of  the  local  coordinates,  it 
is  necessary  to  find  an  expression  which  relates  derivatives  with  respect  to 
the  local  coordinates  to  derivatives  with  respect  to  X_.  This  can  be  done 
using  the  chain  rule  as  follows: 


3X] 

ax2 

H 

H 

3? 

3X1 

ax2 

ax3 

W 

W 

3y 

ax1 

ax2 

CO 

X 

CT> 

az 

az 

az 

ir  > 

9*2  / 


(3.22) 
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(3.13c) 


N3U)  =  \  5  U  +  1) 


In  the  finite  element  approximation  the  displacements  within  each  element 

are  assumed  in  terms  of  the  shape  functions  as  follows: 

3  3  3 

°u  =_^  Ni  (  C)  °u0.  +  y  ( c)  (0£2_£2)i+  Z  1]  Ni (  { °-3"-3) i  (3,14) 


3  3 

au  =  l  N-(c)  Air.  +  i  Ms)  (y  D-  +  z  EJ  A0, 
-  i=i  1  "°i  i=l  1  1  11 


(3.15) 


Figure  4  shows  the  18  degrees  of  freedom  for  the  3  node  beam  element  used 
here.  The  subscript  i  refers  to  the  appropriate  node  and  the  displacements 
are  expressed  in  terms  of  the  nodal  displacements.  Thus,  for  i  =  1,2,3, 

°*o,  -  -Ji.  <v 

*  %  (V 

“ii  ■  *2  (Xo, ) 

s(  -e<‘v  “82.,  °93.) 

Si  ■  E  C8,,.  «e2l,  ”e3l) 

Now  the  displacements  need  to  be  expressed  in  the  form  of  equation  (2.10). 
To  do  this,  define  the  nodal  incremental  displacement  vector  q  as  follows: 


AiT  =  La£oi»  Aiio3’  a!3J 


(3.16) 


Rewriting  equation  (3.15)  in  terms  of  equation  (3.16)  and  using  matrix 
notation. 


*•  '  ,  •  ,  *  ,  *  ^  % 


expressions  in  eq.  (3.90)  are  no  higher  than  4th  order  polynomials,  a 
three-point  numerical  integration  would  give  exact  results.  However,  these 
beam  elements  allow  transverse  shear  deformations  similar  to  a  Timoshenko 
beam.  It  is  well  known  that  when  the  finite  element  approximation  is 
applied  to  such  beams  which  have  a  large  length  to  thickness  ratio  that 
exact  integration  results  in  overly  constrained  equations  which  leads  to 
inaccurate  results  [36],  One  remedy  to  this  situation  is  to  use  reduced 
integration  [29,30].  Therefore  a  two  point  integration  (n  =  2)  was  chosen 
for  eq.  (3.90).  For  n  =  2  the  Gaussian  integration  points  are  aj  =  -1  hfi> 
and  a^  =  1  // 31  whi le  Hj  =  =  1 .  Also,  from  eq.  (3.33) ,  j  |  =  i/2  .  Thus, 

eq.  (3.86)  becomes 


K.  =  |  ^  { T A  (B*)T  C_*  +  \zz  (B*)T  £*  bJ  +  Iyy  (B*)T  c*  (3.91) 


The  initial  stiffness  matrix  and  force  vector  are  calculated  in  a  simliar 
manner.  If  eq.  (3.79)  is  substituted  into  eq.  (3.48)  and  the  integrals  are 
evaluated  as  in  eq.  (3.85),  the  following  results  are  obtained  (again  using 
2-point  numerical  integration): 


1  r 

2  J  =  1 


4  J 


(3.92) 


Finally,  if  eq.  (3.83)  and  (3.77)  are  substituted  into  the  expression  for 
the  initial  stress  force  vector  given  in  Eq.  (3.49)  and  the  integrals  eva¬ 
luated  as  above,  the  result  is  as  follows: 


7  a!a  (lV  °s 


7 

2  j=l 


’sn  +  i  (bJ)t  °su  +  i 

-o  zz  — T  -y  yy 


(Bp'^l^a.  (3.93) 
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CHAPTER  IV 


NUMERICAL  TIME  INTEGRATION  SOLUTION  SCHEME 


4.1  Introduction 

There  are  many  methods  which  can  be  used  for  the  numerical  time 
integration  of  equations  (2.17),  (2.22),  and  (2.75).  Generally  speakiny, 
these  methods  can  be  classified  as  either  explicit  or  implicit  schemes. 
For  any  of  these  schemes  the  time  variations  of  q  and  q  are  discretized 
using  an  m-step,  one-derivative,  linear  multistep  difference  operator 


En+1 


P  P 

ii.  a<  *-■  *  “  i=E., 


b 


(4.1a) 


in+1 


m 

=  i 

i=o 


c .  q  .  + 
l  _n-i 


At 


m 

l 

i=-l 


d .  q 
i  Zn-i 


(4.1b) 


where  qn_i  =  q(tf|_i ) 

Vi  =i(tn-t> 

Vi  *a<Vi> 

At  =  time  increment  =  t  , ,  -  t 

n+l  n 

In  any  particular  application  of  equations  (4.1),  any  of  the  scalar 
constants,  ai ,  b^ ,  c^  or  d-  may  be  zero.  If  both  b_^  and  d  1  are  zero  then 
the  solution  scheme  is  classified  as  explicit.  If  either  of  b_^  or  d  ^  are 
nonzero  then  the  scheme  is  an  implicit  one  and  equations  (4.1)  will 
generally  be  solvable  only  by  an  iterative  procedure  unless  the  differen¬ 
tial  equation  is  linear.  In  this  case  it  is  possible  to  solve  explicitly 


Explicit  schemes  are  computational ly  simplier  than  implicit  schemes 
since  they  don't  involve  the  solution  of  an  implicit  equation.  However, 
with  explicit  methods  it  can  be  shown  [31,32]  that  the  stability  of  the 
solution  is  dependent  upon  the  time  increment  At.  For  linear  systems,  the 
maximum  time  step  for  a  stable  solution  is  related  to  the  highest  natural 
frequency  of  the  system  as  follows: 

At  <  r/»max  (4.2) 

Here  y  is  a  scalar  constant  dependent  upon  the  a^  -  d.  and  w  is  the  hiqhest 
frequency  of  the  system.  For  nonlinear  systems,  the  natural  frequencies  of 
the  system  generally  change  with  time  as  the  configuration  of  the  system 
changes.  Equation  (4.2)  however,  will  still  be  true.  For  certain  "stiff" 
systems  with  widely  varying  natural  frequencies,  the  time  step  necessary 
for  stability  can  be  so  small  that  it  becomes  computationally  time 
consuming. 

Implicit  schemes,  on  the  other  hand,  generally  have  much  less 
stringent  stability  requirements.  In  fact,  some  methods  (e.g.  Newmark's 
method)  are  unconditionally  stable  with  no  time  step  limitations.  The 
accuracy  of  the  solution  then  becomes  the  only  limitation  on  the  time  step 
size. 

For  the  example  problems  discussed  in  this  research,  a  three  node 
beam  element  is  used.  For  small  displacements,  this  is  essentially  a 
Timoshenko  beam.  The  natural  frequencies  for  a  linear  cantilevered  beam 
modelled  by  these  elements  are  shown  in  Tables  1  and  2  as  a  function  of  the 
shear  factor  chosen.  As  is  obvious  from  Tables  1  and  2,  using  these  beam 
elements  will  result  in  a  stiff  system.  The  very  high  frequencies  are  asso¬ 
ciated  with  the  shear  modes  of  the  beam  and  even  through  they  are  only  a 


small  component  of  the  beam  response,  they  control  the  stability  of  the 
numerical  solution.  Even  when  very  small  shear  factors  are  used,  the  high 
frequencies  due  to  shear  will  result  in  unacceptably  small  time  steps  for 
stability.  Thus  even  though  we  are  interested  in  nonlinear  solutions,  the 
linear  solution  clearly  indicates  that  an  implicit  scheme  is  necessary. 

Probably  one  of  the  most  strai ght-forward  and  widely  used  implicit 
schemes  is  the  Newmark  method  or  trapezoidal  rule  [33],  However,  most 
structural  dynamics  problems  involve  second  order  differential  equations 
which  include  only  symmetric  matrices  whereas  the  problem  here  includes  a 
nonsymmetric  gyroscopic  matrix  £.  If  the  Newmark  (or  any  other  strictly 
implicit  method)  is  used  it  would  result  in  having  to  invert  a  nonsymmetric 
matrix  which  involves  a  good  deal  more  computational  time  than  inverting  a 
symmetric  matrix.  In  addition,  many  of  the  terms  in  the  linear  and  angular 
momentum  equations  are  nonlinear  ones  involving  the  displacements.  It 
would  be  preferable  to  treat  these  terms  as  equivalent  "force"  terms  to  be 
moved  to  the  right  hand  side  of  the  equations.  So  what  is  really  needed  is 
a  solution  scheme  which  treats  some  of  the  terms  in  the  equations  expli¬ 
citly  and  some  of  them  implicitly.  Such  a  scheme  is  the  implicit-explicit 
split  operator  method  [25],  In  this  method  terms  such  as  the  £  matrix  are 
treated  explicitly  and  terms  which  control  the  stability  of  the  solution, 
such  as  the  elastic  stiffness  matrix,  are  treated  implicitly. 

4.2  The  Time  Integration  Method 

Consider  again  equations  (2.17),  (2.22),  and  (2.75)  with  all  the  terms 
to  be  treated  implicitly  kept  on  the  left  hand  side  and  all  explicit  terms 
moved  to  the  right  hand  side  of  the  equations. 

H(a)  £  +  P(g)  V0  +  H(q)  £  +  KE(g_)  =  Fq  (R^V^n.^)  -  £U,a)g.  (4.3) 
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P  (q)  q  +  M  V0  +  G(q)  £  =  fT  (%,l0,ntq,q) 

HT(g_)  £  +  GT(£)  V0  +  I_y(q)  £  =  £t  (RQ.V0,£q,q) 


(4.4) 

(4.5) 


Now,  replace  these  three  equations  by  a  single  "general ized"  equation  as 
follows: 


►!g  Sg  *  hi  43g  -  fs 


(4.6) 


where 


(4.7) 
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0 
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(4.9) 


L,- 


£q  ^Sg*  3g^ 

It  (2g»  Se) 
£1t  ^9g»  3g^ 


£  ^g*  Sg*  3 


(4.10 


The  variable  d  is  a  "translational"  variable  the  components  of  which  are 
the  time  integrals  of  V  .  Similarly,  £  is  a  "rotational"  variable  the  com 
ponents  of  which  are  the  time  integrals  of  n.  Thus, 


d 


Assuming  all  the  variables  in  eq.  (4.6)  are  known  at  time  t  = 
equation  at  time  tp+^  =  t  +  At  is  written  as 

Mn+1  *»n+l  ,  ..n+1  .  n+1  rn+l 

_ G  3g  +  £ge  A9g  =  Eg 


tn,  the 


(4.13) 


where 


=  %  (tn+l* 


(4.14) 


=  Hg  (af1) 


(4.15) 


=  Hge  <s2  1 ) 


(4.16) 


-  r  /„n+l  on+1  \ 
'  Ig  *  ^G  ^ 


(4.17) 


using  the  trapezoidal  rule. 


n+1  n  ,  At  / *n  ,  ‘n+1 \ 
3g  ~  9g  +  (3g  +  Hg  ) 


(4.18) 


'n+1  *n  .  At  ,«*n  **n+l . 

Sg  =  3g  +  T  (3g  +  3g  ) 


(4.19) 


Combining  eqs.  (4.18)  and  (4.19)  results  in  the  following: 


_(q"+1  .  q")  __iq"  -%q" 

1 2  G  V  At  JG  — G 


(4.2U) 


Substituting  eq.  (4.20)  into  eq.  (4.13)  yields 
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-G  (it)2  -® 


n,  4  *n  *,n1  .  ,, 
%)  -  IT  %-  MrJ  +  K 


At  2G 


*n+l  Ann+1 

^GE  A^G 


Fn+1 
=  -G 


(4.21) 


The  above  equation  is  the  linearized  form  of  a  nonlinear  equation  and  can 
be  solved  by  Newton-Raphson  iteration.  Letting  subscript  i  represent  the 
i-th  iteration  number,  eq.  (4.21)  becomes 
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Mn+1  r  4  /_n+l  _nx  4  ’_r\  "n,  .  ,.n+l  n+1 


cn+l 
=  £si 


where 


.n+1  n+1  n+1 

%  "  Vl  ’  \ 

Rewriting  eq.  (4.22)  a  bit  results  in  the  following: 


(4.22) 


(4.23) 
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Mn+1  r  4  ,  n+1  n+1  ,  n+1  nx 

•*> 


4  *n  “n-, 

At  2fi  "  Sg3 


,  „n+l  -  n+1  n+K  cn+l 

+  5*,  ‘sw  a®, 1 ' 


(4.24) 


Then,  after  rearranging  a  few  terms,  eq.  (4.24)  becomes 


n+1 i  .„n+l 


[Kqc  +  — Aqf! 
“GEi  (At)2  -Gi  J  ^i 


cn+l  Mn+1  r  4  ,„n+‘ 

F  -  [ - -  (q„ 

"bi  "S’  (At)2  -^i 


4  >  *nni 

--Sg-5g1 


(4.25) 


Using  eq.  (4.20)  and  rewriting  some  of  the  terms  a  bit,  eq.  (4.25)  becomes 


3,  ■  % 


(4.26) 


where  Kt,  =  K^p1  + — ~  m"+1 
_GEi  (At)2"Gi 


(4.27) 


c*  _  cn+l  Mn+1  ••n+1 
fGi  =  £oi  -  -Gi  3g. 


(4.28) 


Thus  by  the  use  of  the  trapezoidal  rule  the  dynamic  equation  of  eq.  (4.13) 
has  been  reduced  to  the  static  equation  of  eq.  (4.26).  The  K£  matrix  is 
called  the  effective  stiffness  matrix  and  is  the  effective  force  vector. 
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The  iteration  starts  with  i  =  1  and  proceeds  until 
required  tolerance,  that  is,  until 


converges  to  within  a 


(4.29) 


where  Tol  is  a  given  tolerance.  For  the  problems  considered  in  this  work, 
Tol=0.01 . 


4.3  Starting  the  Iterations 

The  difficulty  that  now  arises  is  how  to  start  the  iterations.  The 
simplest  way  of  doing  this  is  to  use  the  state  of  the  structure  at  time  t 
as  a  first  approximation  to  the  state  of  the  structure  at  time  tn+1.  Thus, 


(4.30) 

(4.31) 


(4.32) 


This  approximation  is  fine  for  the  mass  and  stiffness  matrices  since  they 
change  slowly  with  time,  but  what  about  the  terms  in  f"+1  which  involve 
and  None  of  these  terms  have  been  treated  implicitly  using  eq. 

(4.18)  as  was  done  for  the  acceleration  and  elastic  stiffness  terms.  The 
reason  for  this,  as  stated  previously,  was  because  many  of  these  terms  are 
nonlinear  or  involve  nonsymmetric  matrices.  In  addition,  the  stability  of 
the  solution  cannot  be  guaranteed  if  these  terms  are  predicted  using  eq. 
(4.30)  [25].  To  ensure  stability,  reference  [25]  indicates  that  when  a 
trapizoidal  rule  is  used  for  the  implicit  part,  the  following  explicit 
"predictors"  for  use  in  the  first  iteration  of  F^ 


must  be  used: 


(4.33) 
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(4.34) 


Tims,  Zg*1  *  £g  (9g+1 •  aS*1 )  <4-35> 

Note  that  eqs.  (4.34)  and  (4.35)  are  obtained  by  setting  *q^  =  0  in  eqs. 
(4.18)  and  (4.19).  Eqs.  (4.30)  -  (4.32)  and  (4.35)  thus  are  used  in  the 
first  iteration  of  eqn.  (4.26).  The  entire  iteration  process  is  summarized 
in  Table  3.  The  solution  scheme  described  in  the  previous  pages  used  the 
regular  Newton-Raphson  iteration.  A  modified  Newton-Raphson  iteration  in 
which  the  matrix  is  kept  constant  throughout  the  iterations  could  also 
be  used.  This  could  result  in  a  computational  savings  since  need  be 
decomposed  only  once;  however,  a  larger  number  of  iterations  will  generally 
be  necessary  than  for  the  reyular  Newton-Raphson  iteration. 

4.4  Spacecraft  Angular  Orientation  Determination 

In  order  to  determine  the  orientation  of  the  spacecraft  at  a  given 
time,  it  is  necessry  to  determine  the  angular  transformation  between  the 
inertial  axes  and  the  rotating  body  axes.  As  shown  in  a  previous  section, 
this  transformation  can  be  represented  in  many  different  ways,  including 
Euler  angles  or  space-three:  1-2-3  angles.  Let  this  set  of  three  orien¬ 
tation  angles  be  represented  by  the  vector  ©  which  is  of  course  a  function 
of  time.  These  orientation  angles  are  related  to  the  body  axis  angular 
velocity  as  follows  [34]: 

0  =  Mg  (9)  £  (4.36) 
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The  MQ  matrix  is  a  3  x  3  matrix  whose  elements  are  functions  of  the  orien¬ 
tation  angles  0.  For  space-three  angles. 
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Euler  angles. 

(4.37) 
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(4.39) 
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(4.40) 


Notethat  for  time  periods  when  cos02  =  0,  space-three  angles  cannot  be  used 
and  Euler  angles  cannot  be  used  if  sine  =  0.  Thus  it  is  necessary  to 
calculate  both  space-three  and  Euler  (or  any  other  description)  so  that  if 
one  formula  fails  at  a  given  time,  the  other  can  be  used. 

It  is  necessary  to  solve  eq.  (4.36)  numerically  and  this  can  be  done 
using  the  following  integration  scheme: 


Vl 


=  ®n  +i7?i0  (0n 


(4.41) 


Both  nn  and  will  have  been  previously  calculated.  Once  the  new  0  has 
been  calculated,  it  becomes  possible  to  determine  the  direction  cosine 
matrix  Cq  which  relates  the  inertial  coordinates  to  the  body  axis  coor¬ 


dinates  as  fol lows: 
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The  elements  of  the  Cg  matrix  for  space-three  angles  are 

gi  ven 

below: 
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(4.43a) 
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• 

Cdj2  =  sin0l  s''n02  cos03  ~  sin03  cos0i 

(4.43b) 

CQ  =  cos©-|  sin©2  cos©3  +  sin©3  sin©^ 

(4.43c) 

- 

• 

CD2,  =  cos02  sin03 

(4.43d) 

C°22  =  S1n0l  sin02  sin03  +  cos03  cos©i 

(4.43e) 

.  • 

C°23  =  COS01  sin02  sin03  ‘  cos03  sin01 

(4.43f ) 

cD31  =  -sin02 

(4.43g) 

**,  •  * 

„• 

1 

cD32  -  sin©!  cos©2 

(4.43h) 

*,t-  * 

cD33  =  cos©1  cos©2 

(4.431 ) 

<r.~  . 

4 

The  elements  of  the  CQ  matrix  for  Euler  angles  are  also 

gi  ven 

as  follows: 

CDll  =  cos  $  cos  ip  -  s  i  n  <j>  cos9  s  i  n  ip 

(4.44a) 

• 

=  -s i n4>  cosij)  -  cos<j>  cose  si n^» 

(4.44b) 

CDn  =  sine  si  nip 

(4.44c) 

„ 

^°21  =  C0S<*>  +  C0S9  C°S ^ 

( 4 . 44d ) 

• 

C^22  =  —si  n  s  i  nip  +  cos  <p  cose  cos  ip 

(4.44e) 

^°23  =  "sin9  cos* 

(4.44f ) 

-  -t 

.  « 

47 

**  ,‘  * 

jmL  n  —  A 

CD31 

=  sin<j>  sine 

(4.44g) 

C°32 

=  cos 4>  sine 

(4.44h) 

C°33 

=  cose 

(4.44i ) 

Once  eq.  (4.41)  has  been  solved  for  0,  then  the  Cq  can  be  calculated  using 
either  Eq.  (4.43)  or  (4.44)  depending  upon  whether  space-three  or  Euler 
angles  have  been  computed  in  eq.  (4.41).  The  CQ  matrix  can  then  be  used  to 
determine  the  orientation  angles  that  weren't  calculated  in  eq.  (6)  [34], 
For  example,  suppose  eq.  (4.41)  is  used  to  solve  for  the  space-three  angles 
and  then  eqs.  (4.43)  are  used  to  calcualte  Since  must  be  the  same 

for  any  orientation  angle  used,  the  Euler  angles  could  be  calculated  by 
solving  for  \p,  e,  and  $  in  eqs.  (4.44).  The  reverse  is,  of  course,  also 
true,  i.e.,  given  CQ,  the  angles  0^  e2,  and  ©3  can  be  solved  for  in  eqs. 
(4.43).  Table  4  summarizes  the  entire  procedure  as  well  as  providing  more 
details  on  solving  for  the  orientation  angles  given  in  the  matrix. 
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CHAPTER  FIVE 


NUMERICAL  RESULTS 

5.1  Gravitational  Forces  and  Moments 

The  method  developed  in  the  previous  sections  is  now  used  to  analyze 
several  example  problems.  The  first  four  examples  involve  the  analysis  of 
structures  which  are  rotating  with  large  enough  angular  velocities  that 
most  of  the  response  occurs  in  a  short  time  period  of  a  minute  or  less. 
Examples  2-4  involve  spacecraft  which  rotate  in  this  fashion.  For  these 
spacecraft,  the  moment  caused  by  the  gravitational  gradient  is  small  and 
can  be  neglected.  Also,  the  spacecraft  will  move  very  little  along  the 
orbital  path  during  the  short  time  periods  which  are  of  interest  in  these 
problems.  The  spacecraft  is  thus  treated  as  though  it  is  in  a  force  free 
environment  (except  for  any  applied  control  forces).  In  example  5,  the 
geometry  of  the  spacecraft  is  such  that  moments  due  to  the  gravitational 
gradient  are  important  and  the  time  periods  of  interest  are  relatively 
long,  up  to  an  hour.  Thus  it  is  necessary  to  obtain  an  expression  for  the 
force  and  moment  due  to  gravity. 

In  Fig.  1,  let  the  inertial  axes  be  located  at  the  center  of  the  earth 
C  and  let  be  the  position  vector  of  the  orbiting  spacecraft.  Neglecting 
any  small  changes  in  R  due  to  deformations, 


where 

yQ  =  Gravitational  constant  of  the  Earth 
M  =  Spacecraft  mass 

The  element  body  force  term  eq.  (2.43)  in  the  principle  of  virtual  work, 
thus  becomes  as  follows: 


“d  =  "  ^d*=  DamPed  natural  frequency 

I33  =  Total  spacecraft  moment  of  inertia 

For  a  given  I^.  the  values  of  the  control  constants  are  determined  by  the 
values  of  w  and  £d  which  are  in  turn  determined  by  the  desired  percent 
overshoot  and  rise  time  of  the  response, 


in2  p 


(- 


(5.23) 

(5.24) 


where  p  =  Percent  overshoot  of  the  •’esponse 

tp  =  Response  rise  time;  the  time  at  which  e  first  equals  eref 

For  p  =  0.01  and  tp  =  5  sec,  the  following  results: 

Sd  =  0.82609 
o)  =  0.902472  rad/sec 
Kt  =  0.814456  I3 
Kw  »  1.49103  I33 

Using  the  spacecraft  parameters  given  in  figure  20  results  in  1 33  =  15,122 

slug-ft2,  Kj  =  12,316  ft-lb/rad,  and  «w  =  22,547  ft-lb-sec.  A  value  of 

0Ref  =  20°  was  chosen  as  the  reference  angle. 

Since  the  spacecraft  is  symmetric,  only  half  of  it  need  be  modeled. 

The  motion  was  studied  using  both  one  and  two  beam  elements  with  a  time 

step  size  of  At  =  0.1  second.  In  addition,  two  values  of  Young's  modulus 

were  used,  E  =  E„  and  E  =  0.5  E  .  The  rotation  angle  for  both  the  rigid 
0  0 


5.5  Example  4  -  Spacecraft  Rotating  to  a  Specified  Orientation 
In  examples  2  and  3  the  motion  of  a  rotating  spacecraft  with 
prescribed  angular  velocity  was  analyzed.  In  this  example  the  more 
realistic  problem  of  a  spacecraft  rotating  to  a  specified  angular  orien¬ 
tation  using  an  applied  control  moment  is  considered.  The  spacecraft  is  the 
simple  symmetric  dipole  pictured  in  figure  20  and  is  originally  at  rest 
until  a  control  moment  that  is  proportional  to  the  angular  change  desired 
is  applied  at  the  center  of  the  rigid  mass.  To  provide  damping,  the 
control  moment  is  also  proportional  to  the  angular  velocity.  The  control 
moment  is  thus  as  follows: 

m  =  -kt  (e-eref)  -Kw  n 

where  =  Control  constants 

6  =  Rotation  angle 

0  f  -  Reference  angle  through  which  the  spacecraft  is  to  be 
rotated 

If  a  sensor  or  some  other  form  of  instrument  for  which  pointing  accuracy  is 
important  is  attached  to  the  spacecraft,  it  becomes  important  to  have  a 
good  estimation  of  how  e(t)  varies  with  time. 

The  rigid  body  solution  for  this  problem  is  the  usual  damped  sinu¬ 
soidal  response  typical  of  underdamped  second-order  control  systems. 


-Cdut  E, 

eRigid  =  eRef  -  e  tcos  wd  1  +  sin  “d  tH  (5*22) 

"  h 

where  w2  =  ^33  =  Circular  frequency 


Cd  =  Kw/2toI 23  =  Damping  constant 


fi,  =  n  sin  y  cos  x  (t-T  ) 

•3  0  0 


(5.20c) 


where  x  =  (— - — )  nn  cos  y  (5.21) 

I  ° 

l22 

For  the  elastic  spacecraft,  the  inertias  of  the  rigid  body  are  chosen  such 
that  1 22  =  I33  >  in  undeformed  structure.  One  beam  element  is  used 
to  represent  each  beam  and  a  time  step  size  of  At  =  0.1  second  was  used  in 
the  numerical  integration. 

The  angular  velocities  of  the  elastic  spacecraft  are  compared  to  the 
rigid  body  case  in  figures  17-19.  The  elastic  response  for  is  much 
closer  to  the  rigid  body  case  than  the  responses  for  the  and  n3  components 
are.  The  three  dimensional  elastic  deformations  are  coupled  with  each 
other  as  well  as  with  the  angular  velocities  causing  a  very  complicated 
response.  The  state  of  the  spacecraft  at  any  given  time  would  be  very  dif¬ 
ferent  from  that  predicted  by  rigid  body  dynamics  alone. 
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6.4  Example  3  -  Three  Dimensional  Precessing  Spacecraft 

In  this  example  the  same  spacecraft  as  in  example  2  is  subjected  to 
out  of  plane  rotations.  As  pictured  in  Fig.  16,  the  spacecraft  is 
"spun-up"  about  a  line  in  the  x^-x^  plane  by  applying  a  moment  during  the 
time  period  0  <_  t  <_  TQ  which  results  in  the  following  angular  velocity: 


9.q  cos  Y  [6(y-)5  -  15(y-)4  +  10  (y-)3]  rad/sec  0  <_  t  <_  T 
o  o  o 


Free  Variable 


t  >  T, 


(t)  tany  0  <_  t  <_  Tc 


Free  Variable  t  >  T, 


0  0  <_  t  <_  T0 

Free  Variable  t  _>  TQ 


(5-1 9a ) 


(5.19b) 


(5.19c) 


where  y  =  Angle  between  the  x-j  axis  and  a 


Values  for  nQ,  Tq  and  y  used  in  this  example  were  as  follows: 

=  0.72552  rad/sec 
T0  =  2  sec 
y  =  30° 

The  rigid  body  response  for  t  >  TQ  in  such  a  situation  would  be  for 
the  angular  velocity  vector  to  precess  about  the  x^  axis  with  a  constant 
precession  rate  A  [35],  For  rigid  bodies  with  I ^  =  1 33  >  Ijj,  the  angular 
velocities  are  as  follows  for  t  >  TQ: 


fil  =  no  cos  Y  =  constant 


(5.20a] 


=  nQ  sin  y  sin  x  (t-TQ) 


(5.20b) 
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The  effect  of  the  control  forces  on  the  rotation  angle  is  illustrated 
in  Fig.  13.  When  a  control  force  is  applied,  the  motion  will  eventually 
damp  out  and  the  average  angular  velocity  will  approach  zero.  As  seen  in 
Fig.  13,  the  larger  the  control  force,  the  smaller  the  rotation  angle. 
Figures  14  and  15  show  the  tip  displacements  and  the  angular  velocity  as  a 
function  of  time.  Note  that  the  amplitudes  of  vibration  decrease  with  time 
as  the  control  forces  are  applied  as  would  be  expected.  One  thing  should 
also  be  noted  in  Fig.  15.  At  t  =  TQ  =  2  sec  the  applied  moment  which  is 
turning  the  spacecraft  is  suddenly  released.  This  sudden  change  or  discon¬ 
tinuity  in  applied  moment  results  in  a  discontinuity  in  the  accelerations 
which  shows  up  as  a  “kink"  in  the  angular  velocity  at  t  =  2  sec.  No  "kink" 
is  seen  in  the  displacements  since  the  velocities  are  all  continuous. 

In  this  example  it  has  been  demonstrated  that  it  is  indeed  necessary 
to  take  into  account  the  deformations  of  the  spacecraft  to  accurately 


Another  way  of  determining  the  effect  of  offset  is  to  see  how  asym¬ 
metrical  the  motion  is.  In  Fig.  11  is  plotted  the  tip  displacement  of  both 
left  and  right  beams  for  IR  =  104  slug-ft2,  mR  =  100  slug  and  Tg  =  15  ft. 

The  displacements  can  be  seen  to  be  almost  the  same  for  each  beam  indi¬ 
cating  that  at  least  for  the  £_(t)  given  in  eq.  (5.16)  the  motion  is  almost 
symmetrical . 

Next,  the  effects  of  these  changes  on  the  rotation  of  the  body  axes 
are  considered.  The  effects  on  the  rotation  angle  of  increasing  the  rigid 
inertia  has  already  been  considered  in  Section  5.2.  The  effects  of 
increased  offset  are  illustrated  in  Fig.  12  which  gives  the  rotation  angle 
as  a  function  of  time  for  rfi  =  0  and  15  ft.  Note  that  for  both  cases  the 
average  angular  velocity  is  the  same  but  that  the  period  is  larger  for  the 
larger  offset. 

5.3.2  The  Effect  of  Control  Forces 

In  this  section  the  effect  of  control  forces  applied  at  the  tips  of  the 
beams  are  considered.  These  forces  are  proportional  to  the  tip  velocity: 

Fc  =  -  Ky  _Vt  (5.17) 

where  Ky  =  Control  constant 
Vj  =  Beam  tip  velocity 

Note  that 

■¥r  =  —  x  — T  +  — T  (5.18) 

£y  =  Position  vector  of  beam  tip 

• 

Uy  =  Elastic  velocity  of  beam  tip 

Three  different  values  of  control  constants  are  considered,  Ky  =  0,  5,  and 
20  Lb-sec/ft.  For  all  three  cases,  the  rigid  mass  and  inertia  were  mR  = 

500  slugs  and  IR  =  5( 1 0) 4  slug-ft2  with  an  offset  of  rfi  =  5  ft.  Once 
again,  a  time  step  of  0.1  sec  was  used  in  the  numerical  integration. 


5.3.1  Parametric  Studies 

For  all  the  problems  in  this  part,  the  beam  properties  were  kept 
constant  and  only  the  properties  of  the  central  rigid  mass  were  varied. 

The  way  these  changes  affect  the  frequency  of  vibration  of  the  beams  can 
thus  be  determined. 

The  equations  of  motion  were  solved  using  a  time  step  of  At  =  0.1  sec. 
and  one  element  to  represent  each  beam.  The  resulting  beam  periods  of 
vibration  are  summarized  in  table  5.  Note  that  the  period  increases  as  the 
size  of  the  rigid  mass  increases.  As  the  rigid  inertia  grows  smaller  in 
comparison  to  the  inertia  of  the  beam,  the  problem  approaches  that  of  a 
free-free  beam  vibrating  in  its  second  mode  shape.  For  a  200  ft.  length 
beam,  the  linear  second  mode  frequency  results  in  a  period  of  2.2  seconds. 
As  the  rigid  inertia  grows  very  large  in  comparison  to  the  beam  inertia, 
the  problem  approaches  that  of  a  cantilever  beam  attached  to  a  very  large 
rigid  mass.  The  linear  period  for  a  100  ft.  cantilever  is  6.4  seconds.  In 
Table  5,  the  periods  will  vary  between  these  two  extremes. 

Now,  what  about  the  effect  of  the  offset?  As  can  be  seen  in  Table  5, 
for  a  given  rigid  inertia  as  the  offset  increases  so  also  does  the  period. 
However,  notice  that  the  change  in  period  is  much  less  as  the  rigid  inertia 
increases.  Also,  to  determine  the  effects  of  varying  the  rigid  mass  while 
keeping  the  rigid  rotary  inertia  constant,  the  following  was  done.  For 
IR  =  104  slug-ft2  and  Tq  =  15  ft,  rigid  masses  of  mR  =  100,  500,  and  1000 
slugs  were  used.  The  resulting  periods  were  all  the  same  value  of  T  =  3.3 
sec.  indicating  that  mR  has  little  effect  on  the  period.  What  was  affected 
by  mR  was  the  translation  of  the  body  axes.  These,  however,  were  very 
small  and  never  exceeded  one  foot  and  thus  had  very  little  effect  on  the 
rest  of  the  motion. 
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5.3  Example  2  -  Rotating  Spacecraft  with  an  Offset  Center  of  Mass  and 

Applied  Control  Forces 

In  this  example  a  simple  spacecraft  consisting  of  two  slender  aluminum 
beams  attached  to  a  central  rigid  mass  is  analyzed  (Fig.  10).  The  rigid 
mass  has  a  mass  moment  of  inertia  of  IR,  a  mass  of  mR,  and  a  width  of  2r^. 

A  set  of  body  axes,  x^ ,  x^  and  x^  are  fixed  in  the  spacecraft  and  rotate 
with  it.  The  centerlines  of  the  beams  are  originally  aligned  along  the 
Xj  axis  and  the  center  of  mass  of  the  spacecraft  is  offset  a  distance 
Tg  along  the  x^  axis. 

The  spacecraft  is  originally  at  rest  but  at  time  t  =  0  a  moment  is 
applied  to  the  spacecraft  for  TQ  seconds  such  that  a  given  angular  velocity 
is  achieved.  This  angular  velocity  was  chosen  so  that  the  rigid  mass  would 
rotate  through  a  given  angle  in  the  x^  -  plane  and  then  come  smoothly  to 
rest  at  t  =  TQ.  After  Tq  seconds,  the  applied  moment  is  released  and  the 
angular  velocity  becomes  a  free  variable.  The  particular  form  chosen  for 
the  angular  velocity  was  as  follows: 

f 

-f  M  +  sin  *  (-r-  -  A)  ]  rad/sec  0  <_  t  <_  T 

J  o  c 

n( t)  =  )  (5.16) 

I  Free  Variable  t  >  TQ 

For  all  the  problems  solved  in  this  example,  nQ  =  tt/20  rad/sec  and  Tq  =  2 
seconds.  Note  that  if  rfi  =  0  the  problem  is  symmetric  and  only  one  half  of 
the  spacecraft  need  be  modeled.  However,  if  rQ  +0,  then  in  general  the  mo¬ 
tion  will  not  be  symmetric  and  the  entire  spacecraft  will  need  to  be  modeled. 

The  results  of  this  example  are  divided  into  two  parts.  In  the  first 
part,  a  parametric  study  is  done  in  which  the  effects  of  varying  the  rigid 
mass  and  inertia  as  well  as  the  center  of  mass  offset  are  studied.  In  the 
second  part,  the  effect  of  applied  control  forces  is  examined.  For  both 
cases,  eq.  (5.16)  is  used  during  the  "spin-up"  period  of  motion. 
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will  be 


Also  note  that  the  larger  IR  is  in  comparison  to  Ig,  the  closer  nAvg 

to  the  rigid  body  angular  velocity  nQ. 

Now  eq.  (5.15)  can  be  used  to  explain  the  behavior  of  the  beam  in 

figures  8  and  9.  At  t  =  1.0  sec  the  beam  is  moving  away  from  the  body  axis 

and  thus  Hn  <0.  At  t  =  3.3  sec.  the  beam  has  reached  its  maximum  displa- 
uo 

cement  relative  to  the  body  axis  and  thus  Hn  =0.  At  t  =  7.0  sec.  the 

uo 

beam  is  leading  the  rigid  body  motion  and  Hn  >  0.  Thus  the  average  anyu- 

uo 

lar  velocities  pictured  in  figure  8  are  as  predicted  by  eq.  (5.15). 

Similar  results  are  obtained  in  figure  9  except  that  with  a  larger  inertia 
ratio,  all  the  averge  angular  velocities  are  closer  to  as  predicted 
again  by  eq.  (5.15) . 
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where 


IB  =  Mass  moment  of  inertia  of  the  beam  about  the  shaft 
centerline 


Also,  let 


Hr,  = 


HTi 


(5.11) 


The  value  of  Ig  will  change  only  slightly  as  the  beam  vibrates  and  for  the 
purposes  of  this  discussion  it  is  sufficient  to  consider  it  to  be  a 
constant.  The  conservation  of  angular  momentum  requires  that 


(IR  *  V  %  +  Ho0  -  Ur  *  IB>  0+  Hd 

"here  hd0  -  (hd),.t> 

Solving  for  n. 


(5.12) 


'MOM 


t  =  v  (T^'  <\  -  hd> 


(5.13) 


The  angular  velocity  will  be  constantly  changing,  but  what  will  its  average 
value  be?  Taking  a  time  average  of  eq.  (5.13), 


!Avg  %  +  (IR  +  iB>  -  (^Avgl 


(5.14) 


Since  the  beam  vibrates  about  the  equilibrium  position,  ( hd) Avg  ^  be 
small  and  can  be  negleccted.  Thus, 


So  if 


fiAvg  % 

-r- )  HD 

[B  0 

Hn  <  0 
u0 

then 

nAvy 

Hg  >0 
u0 

then 

nAvg 

o 

II 

o 

o 

31 

then 

nAvg 

(5.15) 
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The  results  for  case  2  are  presented  in  figures  8  and  9,  which  show  the 
variation  of  the  rotation  angle  of  the  body  axes  with  time.  In  figure  8, 
IR/Ig  =  0.5  and  in  figure  9  IR/Ig  =  1.0.  For  both  cases,  TMQM  =  1.0, 3. 3, 
and  7.0  seconds.  In  both  figures,  the  "average"  angular  velocity  is 
greater  than  for  Tu-u  =  7  and  less  than  n  for  =  1.0  while  for 
Tmom  =  3.3  it  is  about  the  same  as  Note  also  that  all  the  angular 
velocities  are  closer  to  for  the  larger  inertia  ratio  in  figure  9  than 
in  figure  8. 

In  order  to  understand  the  behavior  of  the  beam  in  figures  8  and  9, 
consider  what  happens  when  the  moment  support  is  released  at  t  =  TMqm  and 
the  angular  velocity  becomes  a  free  variable.  For  t  >  TMqM,  there  is  no 
longer  an  applied  moment  turning  the  shaft  and  the  angular  momentum  will 
remain  constant.  The  angular  momentum  is  determined  as  follows: 

p(  £  x  dv  (5.6) 

v 


But 


dr 

-ft  =  £x  I  +  iL 


(5.7) 


Substituting  eq.  (5.7)  as  well  as  eq.  (2.14)  into  eq.  (5.6)  and 
rearrangi ng. 


—2 

P  r* 


dv)  Q  + 


p  £  N_*  dv)i^ 


V  V 

Using  eqs.  (2.29)  and  (2.30)  results  in 
^  =  l  U  +  H_T  £ 

For  this  particular  problem,  J^,  I_  and  j^are  scalars  and 

I  -  Ir  +  h 


(5.8) 


(5.9) 


(5.10) 
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Case  2  -  n(t)  free  for  t  >  TV 


n(t)  = 


f6  (T“)5  -  15  (f-)4  +  10  rad/sec  0  1  t  1  T 

0  0  o  u 


nQ  rad/sec 
Unknown  variable 


Toi*i  TM0M 

1 1  thoh 


For  both  cases,  values  of  =  it/10  rad/sec  and  Tq  =  1  sec  were  used.  The 
equations  of  motion  were  integrated  using  a  time  step  of  At  =  0.1  sec  and 
one  beam  element  was  used. 

The  results  for  case  1  are  shown  in  figures  6  and  7.  Fiyure  6  shows 
the  angular  displacement  of  the  tip  of  the  beam  as  it  rotates.  This  is  the 
angle  between  the  inertial  axis  and  a  line  drawn  from  the  origin  of  the 
body  axis  to  the  tip  of  the  beam  at  time  t.  The  body  axes  and  inertial 
axes  coincide  at  t  =  0.  In  figure  7  is  pictured  the  orientation  of  the 
deformed  beam  at  different  times.  The  dashed  lines  represent  the  positions 
of  the  x1  body  axis  at  the  given  times  and  coincides  with  the  motion  of  a  rigid 
beam.  The  solid  lines  are  the  shapes  of  the  deformed  beam  and  the  arrows 
represent  the  displacement  of  the  beam  tip. 

The  linear  natural  frequency  for  this  beam  as  given  in  Table  1  is  wj  = 
0.547  rad/sec  which  corresponds  to  a  period  of  11.5  sec.  Note  that  a>^> 
and  thus  it  would  be  expected  that  the  beam  would  complete  about  1.7 
vibration  cycles  about  the  rigid  body  position  during  the  20.5  sec.  it 
takes  for  the  shaft  to  complete  one  rotation.  As  can  be  seen  in  figures  6 
and  7,  this  is  exactly  what  the  beam  does.  For  the  first  6  seconds  it  lags 
behind  the  rigid  body  motion,  then  for  the  next  6  seconds  it  leads  the 
rigid  body  motion  completing  a  cycle  in  12  seconds.  Thus  the  beam  vibrates 
back  and  forth  about  the  rigid  body  motion. 


5.2  Example  1  -  A  Beam  Rotating  About  a  Fixed  Axis 


In  this  example  a  thin  uniform  aluminum  beam  is  fixed  to  a  shaft  which 
is  rotating  about  its  centerline  with  an  angular  velocity  n.  The  shaft  has 
a  mass  moment  of  inertia  of  IR  about  its  centerline  which  is  fixed  so  that 
no  translations  occur.  The  radius  of  the  shaft  is  negligible  compared  to 
the  length  of  the  beam.  A  set  of  body  axes  x-j ,  x^  and  x3  are  fixed  in  the 
shaft  and  rotate  with  it.  At  time  t  =  0  the  beam  and  shaft  are  at  rest  and 
the  beam  is  aligned  along  the  x-j  axis.  The  situation  is  illustrated  in 
figure  5. 

Two  different  cases  are  considered.  In  case  1  the  angular  velocity  n 
is  a  specified  function  of  time  for  all  time.  The  only  unknowns  are  the 

nodal  displacements  and  velocities  which  can  be  found  using  eq.  (2.75)  with 

•  • 

n_  and  given  and  J/Q  =  ^  =  _0.  For  case  2,  the  angular  velocity  is  spe¬ 
cified  only  for  t  <_  TMQM  after  which  it  becomes  a  free  variable.  The  time 
period  0  _<  t  _<  TMQM  is  the  "spin  up"  period  during  which  the  beam  is  acce¬ 
lerated  from  rest  to  a  given  angular  velocity.  For  t  >  TMqM,  the  angular 
velocity  becomes  an  unknown  variable  and  eq.  (2.28)  must  be  used  to  solve 
for  it.  For  each  case,  the  angular  velocity  during  the  spin  up  period  is  a 
smooth  polynomial.  The  two  cases  can  be  summarized  as  follows: 

Case  1  -  a{ t)  given  for  all  time 
r 

%  [6  (y- )  -  1 5  (y-)  4  +  10  (y— )  3]  rad/sec  0  <_  t  <_  TQ 

0  0  0 


n(t)  = 


rad/sec 


t  >  Tn 


£b  •  -(  /  P  N*1  dv)  — ^  (5.2) 

K  K 

e 

or 

(5.3) 


The  moment  due  the  gravitational  gradient  can  be  shown  [34]  to  be  as 
follows: 


—  5c  rG 
0 


P  (eJ  r) (R0  £)  dv 


where 


p  r  dv 


Note  that  r„  is  the  instantaneous  location  of  the  center  of  mass 


body  case  and  the  elastic  cases  are  pictured  in  figure  21.  Note  that  for  E 
=  Eq,  the  response  time  is  still  about  5  seconds  but  that  because  of  the 
elasticity  of  the  beams,  the  peak  overshoot  is  near  10%.  Also,  it  takes 
about  12  seconds  before  the  orientation  does  not  vary  more  than  1%  from  20°. 
The  more  flexible  beam  with  E  =  |  EQ  still  has  a  peak  overshoot  of  10%  but 
the  rise  time  is  increased  to  5.8  seconds  and  it  takes  even  longer  for  the 


rotation  to  settle  down  to  within  1%  of  0Ref.  Figure  22  shows  the  actual 
orientation  of  the  spacecraft  for  E  =  EQ  at  t  =  5  seconds.  The  beam  tip 
lateral  displacement  are  plotted  in  Figure  23.  Note  that  both  figures 


illustrate  the  large  deformations  which  occur. 
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5.6  Example  5  -  Gravity  Gradient  Stabilization  of  a  Spacecraft  in  a 
Circular  Earth  Orbit 

The  previous  example  examined  spacecraft  with  fairly  high  rates  of 
rotation  during  relatively  short  time  periods.  This  example  considers  a 
slowly  turning  spacecraft  over  a  time  period  which  is  a  significant  percen¬ 
tage  of  the  orbital  period.  Consider  the  spacecraft  pictured  in  figures  24 
and  25  which  consists  of  a  single  long  slender  beam  attached  to  a  rigid 
mass.  The  spacecraft  is  in  a  200  mile  altitude  circular  Earth  orbit  and  at 
time  t  =  0  the  beam  axis  is  at  an  angle  8  =  45°  to  the  orbital  radius  vec- 
tor  Rq.  The  spacecraft  is  rotating  with  an  initial  angular  velocity 
“o  =i%/Ro  =  1  -1512(10-3)  rad/sec.  and  thus  8  =  0.  The  gravitational 
gradient  causes  a  moment  (eq.  (5.4))  which  tends  to  align  the  beam  axis 
with  the  orbital  radius  vector.  A  control  moment  Mc  is  applied  at  point  0 
in  the  rigid  mass  as  follows: 


”c  ’  -Le  6 


where  L.  =  Control  constant 

P 


(5.25) 


In  order  to  determine  the  values  of  Lg  to  use,  consider  the  rigid  body 
problem  with  smal 1  8. 


8  +  2yp8  +  p  8  =  0 


(5.26) 


With  y  =  1,  the  system  becomes  critically  damped,  thus  for  an  underdamped 
system,  choose 

Le  <  2pl33  (5.29) 

The  spacecraft  is  in  a  200  mile  high  orbit  and  thus  the  orbital  period  is 
5458  sec.  Two  different  sizes  of  rigid  mass  are  considered: 

Case  1  Case  2 

mR  =  5000  slugs  mR  =  1500  slugs 

IRl  =  104  ft-lb-sec2  IRl  =  104  ft-lb-sec2 

Ir2  =  108  ft-lb-sec2  IR  =  2.43  ( 1 06)  ft-lb-sec2 

JR3  =  3(1q8)  ^-lb-sec2  IR  =  3(  1 06)  ft-lb-sec2 

The  size  and  material  properties  of  the  beam  are  identical  for  each  case. 
The  problems  were  tested  using  both  a  one  and  two  element  model  with  a  time 
step  of  At  =  10  seconds. 

Figures  26  and  27  show  the  results  for  case  1  for  both  free  vibration 
(no  control  forces)  and  a  control  moment  with  Lg  =  4 .5( 1 04)  ft-lb-sec. 

Note  that  even  though  fig.  25  indicates  that  large  elastic  displacements  in 
the  beam  are  occuring,  the  motion  of  the  rigid  mass  is  almost  the  sames  as 
that  for  a  rigid  body.  Thus,  for  such  a  large  rigid  mass,  the  elastic 
displacements  of  the  beam  have  little  effect  on  the  variation  of  the  angle 


The  results  for  case  2  are  pictured  in  figs.  28  and  29.  For  no 
control  forces,  the  elastic  deformations  in  this  case  do  have  an  effect  on 
B,  even  though  it  generally  follows  the  rigid  body  case.  With  the  control 
force  added,  however,  the  motion  for  e  is  nearly  indentical  to  the  rigid 
body  case. 
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CHAPTER  VI 
CONCLUSIONS 

The  numerical  results  obtained  in  Chapter  V  indicate  that  the  method 
developed  here  is  capable  of  determining  the  time  response  of  unrestrained 
flexible  structures  which  are  undergoing  large  elastic  deformations  coupled 
with  gross  nonsteady  rigid  body  translational  and  rotational  motions  with 
respect  to  an  inertial  reference.  The  use  of  an  implicit-explicit  split 
operator  numerical  integration  scheme  has  resulted  in  stable  solutions  for 
all  the  problems  tested.  In  addition,  the  example  problems  indicate  that 
the  method  is  capable  of  analyzing  problems  which  include  the  effect  of 
control  forces.  Although  only  beam  elements  have  been  used  in  this  work, 
the  equations  in  Chapter  II  are  quite  general  and  will  apply  for  more 
complicated  elements  such  as  plates  and  shells.  Also,  the  method  can  be 
used  to  solve  problems  which  include  more  complicated  types  of  motions  such 
as  spacecraft  deployments  involving  rotations  and  relative  velocities  bet¬ 
ween  different  spacecraft  parts. 
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Table  1 


Natural  Frequencies 

of  a  Linear  Cantilever  Beam  Modelled  by  One 
Timoshenko  Element 

L/t  =  100 

3-Node 

.-•Vv 

• 

Shear  Factor  8 

u-j  <^2 

co_ 

3 

“4 

(Rad/Sec) 

• 

1.0 

0.5470  4.502 

8417 

11, 

287 

0.5 

0.5470  4.499 

5956 

7, 

981 

0.1 

0.5470  4.478 

2677 

3, 

570 

r 

.  • 

0.01 

0.5448  4.259 

891 

1, 

132 

m-  *<■»«* 

Table  2 

Natural  Frequencies 

of  a  Linear  Cantilever  Beam  Modelled  by 
Timoshenko  Elements 

L/t  =  100 

Two  3 

-Node 

Shear 

;•  /■- 

Factor  w, 

8 

w2  w3  “4  “5 

% 

“7 

w8 

(Rad/Sec) 

• 

1.0  0.5222 

3.533  11.88  49.85  5271 

9900 

11,287 

11,294 

0.5  0.5528 

3.549  11.87  49.39  3763 

7007 

7,981 

7,988 

0.1  0.5311 

3.524  11.74  46.10  1808 

3157 

3,570 

3,582 

- 

0.01 


0.5312 


3.411  10.62  29.94  905 


1073  1,132  1,173 


TABLE  3 


Summary  of  Newton-Raphson  Iteration  Scheme 

—  -  <• 

Predict  quantities  for  first  iteration,  i=l 

• 

nn  +  l 

n 

= 

• 

*n+l 

*n 

=  Sg 

‘  '  ' 

••n+1 

4  >  "n 

m 

*, 

At  -Sg  "  3g 

W 

as;1 

„n  .  *n  .  (At) 2  **n 

9g  +  At  Sg  ~  •  Sg 

*n+1 

*n  .  At  “n 

q-G, 

=  Sg  +  T  Sg 

■  .*  '  " 

Form  Modified  Stiffness  Matrix  and  Force  Vector 

•  ] 

II 

sf 

li 

c  ^n+l  “Tn+lv  M  /rtn+lx**n+l 

Ig  ^3g1  *  Sg1  ^  "  !!g  ^Sg1  )  Sci 

Solve  for 

incremental  displacements 

• 

.  n+1 

*1  * 

^i 

=  Kp .  Ff. 

-G-l  — Gi 

■  • 

Update  or 

correct  the  displacements,  velocities  and  accelerations 

• 

n+1 

n+1  .  ,  n+1 

Qr 

-ui+l 

=  Sg1  l  Sg1 

'•V*V  * ^ 

p 

••n+1 

3si+l 

4  ,  n+1  n>  4  *n  **n 

(At)2  <3S1  38  4t  %  36 
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TABLE  3 


If  yes,  go  to  step  6 

If  no,  i  =  i  +  1  and  yo  to  step  2 

6.  Update  angular  orientation  (see  Table  4) 

7‘  *0+1  ■  *0+1  *  4* 

*n+l  —  *max,  go  to  step  1 

lf  *o+l  >  W  stop 
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TABLE  4 


SPACECRAFT  ANGULAR  ORIENTATION  DETERMINATION 


Update  the  space-three  angles 
If  0  =  90°,  Go  to  step  4 
Otherwise,  use  eq.  (4.39)  and  (4.40) 


Vl  ■  +  4  Me  <2h  +  Vi  I 


Calculate  CQ  using  eqs.  (4.45) 


Calculate  Euler  Angles 

a)  If  |C33|  =  1 ,  go  to  step  3c 


otherwise 


.-1 


e  =  cos  (C33)  0  <  e  <v 


b)  Let  a  =  cos”1  (-  0  <_  a  <_  it 


*  -  cos_1  (s^e) 


0  <  6  <  ir 


Then 


*  = 


a  if  C  -j  3  >  Q 


2w  -a  if  C13  <  0 


and 


8  if  C3l  1  0 


2ir  -8  if  C31  <  0 


stop 


c) 


r  0  if  c,«  =  1 


TABLE  4 


<t>  =  0 

a  =  COS-1  (C22)  0  <_  a  <_  it 


^  f  aifC21i° 

2 tt  -a  if  C2 j  <  0 

stop 

Update  Euler  Angles 
Use  eqs.  (4.41)  and  (4.42) 

-n+l  *  2n  *4*0  (Sn  *  fihtl  > 


Calculate  Cq  using  eqs.  (4.46) 

Calculate  Space-Three  Angles 

a)  If  |c^|  =  1 ,  go  to  step  6c 
Otherwise, 

02  *  si n  1  (^31 )  ”  7  —  ®2  —  7 

b)  Let  a  =  sin  (cos'e"1  ~7<-a<.J 


6  =  sin"1  ( 


'21 


cos  ©„ 


IT  IT 

7  -  0  -7 


Then, 


0 


1 


J  01  i  f  >_  0 
ir-a  if  Cjj  >  0 


9 


3 


Stop 


e  if  cn  0 

TT- 3  if  C1 !  <  0 
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TABLE  5 


Period  of  Vibration  of  the  Beams  for  Various  Rigid  Body  Inertias 


offset  rg  (Ft) 


0 

_ 

5 

15 

IR  x  (104),  mR 

Period  of  Vibration  (sec) 

1.0,  100 

2.7 

2.8 

3.3 

2.5,  250 

3.3 

3.4 

3.9 

5.0,  500 

4.1 

4.2 

4.5 

20,  2000 

6.3 

6.4 

6.4 
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Rotate  ij»  about  x3  Rotate  e  about  £  Rotate  about  5 

(b)  Space  -  Three:  1-2-3  Angles 


Figure  6 

Angle  of  Beam  Tip  for  Rotating  Beam 


20.0  r  Fioure  11 

Bearn  Tip  Displacements  for  an  Offset  of  15  Ft. 


Rotation  Annie  b  -  Case 


Figure  25 

Spacecraft  Geometry 


10.0 


(}j)  }u9iu9DeidsiQ  dii 


Beam  Tip  Displacements 


Figure  20 

Example  4  -  Spacecraft  Rotating  to  a  Specified  Orientation 


Figure  16 

Example  3  -  Three  Dimensional  Precessing  Spacecraft 


Rigid  Mass  Properties: 


Ip  =  6.1(1 0^ )  slug-ft 
K11 

Ip  =  5.0(1 04)  slug-ft2 
K22 

Ip  =  4.829175(1 04 )  slug-ft 
K33 


Beam  Properties: 


L  =  100  ft. 

A  =  0.0654  ft2 

1=1=  0.008185  ft* 
yy  zz 

p  =  5,22  slug/ft 
1 . 44 ( 1 08 )  lb/ft2 
5 . 54 ( 1 0 7 )  lb/ft2 


mR  =  500  slugs 


Figure  28 


APPENDIX  A  -  ORIENTATION  MATRICES  AND  VECTORS 


Let  a^,  a^,  and  a^  be  a  set  of  mutually  orthogonal  axes  fixed  in  the 
undeformed  beam,  centered  at  the  centroid  of  the  cross-section  with 
a^  tangent  to  the  longitudinal  axis  of  the  beam.  Also,  °e1  °02,  and 
°e3  represent  the  space-three  1-2-3  orientation  angles  of  the  beam.  Let 

f"  h 

d^j  and  e^.  be  the  ijl  element  of  the  matrices  D  and  £  respectively.  Th 
these  elements  are  as  follows  for  k  =  1,2,3. 


=  (sin°e2  cos°0j  cos°e3  +  sin°0j  sin°03)  a^k  +  (sin°e2  sin' 
cos°e1  -  sin°e1  cos°03)  a2k  +  cos°01  cos°02  a3k 

=  s i n °e ^  cos°02  cos°03  alk  +  sin°0j  sin°03  cos°02  a2k 
sin°02  sin°0^  a3k 


in°02  sin°01  a3k 

dk3  =  ~(sin°01  sin°e2  sin°03  +  cos°01  cos°03)  a]k  +  (sin°01 
cos°03  -  sin°03  cos°0^ )  a2k 


1  °0g  -  sin°03  cos°0| )  a2k 

in°0^  sin°02  cos°03  +  sin°03  cos°0j)  a^k  -  (sin°0 


ekl  =  (■s'’n°9i  sin°02  cos°03  +  sin°03  cos°0j)  a^k  -  (; 
sin°03  +  cos°01  cos°03)  a2k  -  sin°0^  cos°02  a3k 

o  _  _ oa  _ o.  _  _ _ o„  _._o.  _ 


=  cos°0^  cos°02  cos°e3  a-|k 
sin°02  cos°0^  a3k 


+  cos°0^  sin°03  cos°02  a2(( 


ek3  =  (-cos°9]  sin°02  sin°03  +  sir 
cos°03  +  sin°01  sin°03)  a2k 


in°0^  cos°03)  alk 


Also,  the  change  in  directional  axes  is  as  follows: 


For  k  =  1,2,3 
(°a2  -  a2>k 


sin°02  cos°03  -  cos°e^  sin°e3)  a^  + 

^  sin°02  sin°03  +  cos°01  cos°03  -  1 )  a? 

rnc°fl_  a 


(°a3  “  a3)k 


sin°0|  cos°02  a^k 

(cos0©^  sin°02  cos°03  +  sin°0|  sin°03) 
(cos°0^  sin°02  sin°03  -  sin°0^  coso03^ 

(cos°01  cos°02  -  1)  a3k 


alk  + 
a2k  + 


follows: 


and  (A8)  are 

related  to  the  D_  and  E  matrices  as 

-  a2>k  '  d1k 

i,k  =  1,2,3 

( A9) 

'  a3^k  °  e1k 

i,k  =  1,2,3 

(A10) 

Finally,  the  transformation  matrix  T  between  local  and  global  stresses 
and  strains  is  as  follows: 


af, 

4 

2 

al  3 

2a 11 3 1 2 

2al lal 3 

2al 2al 3 

2aUa21 

2a12a22 

2al 3a23 

alla22+a12a21 

alla23+a13a21 

a12a23+a13 

2al la31 

2al 2a32 

2al 3a33 

al  1  a32+al 2a31 

alla33+a13a31 

al 2a33+al 3 

APPENDIX  B 


ELEMENTS  OF  MATRICES  OF  WHICH  THE  STIFFNESS  MATRIX  IS  COMPOSED 

The  Dk  Matrix 

a - - 

Dk  =  Dq( c)  +  y  D*(c)  +  Z  d£u)  k  =  1,2,3 
k  k  k 

The  oo,  D^,  and  Dz  are  partitioned  into  three  3x6  submatrices  as  follows: 
-o  ■  HsJ),  (fi*)2  (fio ) 3 1 

“y  *  [(Dy),  (Oy)2  (oJ)3] 

-  t(a^),  (&z>2  <fi,>3] 

Now,  let  (drs)i  and  (eps)^  be  the  elements  of  the  Di  ad  E_^  matrices 
respectively  (see  Appendix  A)  and  Ni  be  the  ith  shape  function,  then  for 


i,k  =  1,2,3 


3Ni 

0 

0 

0 

0 

0 

m  *  +• 

_  • 

(“Jli  * 

0 

0 

0 

Ni<dn)i 

Ni (d 1 2 ) i 

NJ<dl3>1 

(64) 

0 

0 

0 

¥en>i 

Ni(ei2)i 

Vel3>1 

• 

0 

3N. 

1 

~n 

0 

0 

0 

0 

-  • 

<2o>i  ' 

0 

0 

0 

N,(d  )i 

1  21  1 

Md  ). 
22 

N  (d  , 

23  i 

(65) 

0 

0 

0 

Ni(e2i)i 

Ni ( e22 )  i 

Ni ^  e23  ^ i 

(61) 

(62) 

(63) 
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where 


hrt  = 
Ok 

3 

z 

i*l 

3N 

i  ,o  x 

3E  "0k  ' 

(B12) 

1 

"•  . 

II 

£ 

3 

z 

i=l 

3N. 

__ I  (V  _ 

'  a2k 

3£  C 

a2k^i 

(B13) 

• '  •’  _  " 

h  = 

ozk 

3 

z 

i=l 

3N. 

__ L  /°a«  . 
3C  (  3k 

a3k>1 

(B14) 
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• 

hyk  = 

3 

z 

i=l 

Ni(0  (°a*k 

"  a2k^i 

(B15) 

• 

hzk  = 

3 

z 

i=l 

(°a'k 

“  a3k^i 

(B16) 

Let  the  inverse  of  the  Jacobian  matrix  be  represented  as  follows: 


l 

1 

L13 

•  '  -.**  .  ‘ 

L11 

L1 2 

-• 

i'1  - 

L21 

L22 

L23 

(B17) 

!  *  • 

L31 

L32 

L33 

Note  that  J_1 

will  be  a  constant.  Thus  the 

3° 

derivatives  {— g^-}  can  be  written 

i 

using  eq.  (55) 

as  follows: 

3°uk 

K 

axj 

^kj  +  y  y*kj  +  z  lkj 

k.j  =  1,2,3  (B18) 

• 

.  •  -  ■ 

’  '  o*. 

°*kj 

=  Ljl  hok  +  Lj2  hyk  +  Lj3  hzk 

( B 1 9 ) 

V  \/v 

_• 

„ 

y^kj 

=  LJl  Vk 

(B20) 

**kj 

=  LJl  ho^k 

(B21) 

_• 

Finally,  by  substituting  the  derivatives  in 

eq.  (B18)  into  the  expressions 

no 

„• 

v 

1 

.  ~  -1 

/  V  “  ►  •  '  •  *  -  *  *  "  •  •  •  *  .  •  *  .  •  K  '  ,  '  ' 

APPENDIX  C 


DETAILS  OF  MATRICES  AND  VECTORS  USED  IN  THE  MOMENTUM  EQUATIONS 


LINEAR  MOMENTUM  EQUATION 
The  Element  £  Matrix  (18x3) 


(Cl) 


Substituting  eqs.  (3.18)  and  (3.13)  into  eq.  (Cl)  and  carrying  out  the  volume 
integral  results  in  the  following  18  x  3  matrix  for  P: 

fe  =  "ip  U3  0  4I3  0  I3  0]  (C2) 

where  i  =  Beam  element  length 

I3  -  3  x  3  Indentity  Matrix 
0  =  3  x  3  Null  Matrix 


The  £  Matrix  (3X3) 


(C3) 


pIo 


dv  + 


p  °u  dv  + 


p  Au  dv 


(C4) 


Now, 


p  £ 0  dv  *  Matrix  of  Center  of  Mass 
”  Position  at  t  =  0 


(C5) 


In  many  cases  GQ  =  0.  Substituting  eqs.  (3.14),  (3.17),  (3.18) -( 3.21)  into 
eq.  (C4)  and  integrating  over  the  volume  results  in  the  following: 


Jt 
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where  M  .=  Total  Mass 


Hg  =  irf  p(°“  +  A^dv  =  mIT  (° 


1  +  AgJ 


(C7 ) 


v 

Note  that  is  the  change  in  the  center  of  mass  from  the  initial  unde¬ 
formed  state.  Eq.  (C6)  can  be  simplified  to  eq.  (C8).  Note  also  that  the 
use  of  the  P  matrix  is  after  its  assembly. 


G 


T 


where 


rr  =  r„  +  u~ 
— b  —bo  — b 


(C8) 

(C9) 

(CIO) 


The  f_c  Force  Vector  (3x1) 


p  r  dv  = 


p {£  +  °Au_  +  Au)  dv 


(Cll) 


Using  eq.  (C7)  and  (CIO),  eq.  (Cll)  becomes 


Ic  =  M  J?  IQ 


(Cl  2) 


ANGULAR  MOMENTUM  EQUATION 
The  Matrix  (3x3) 


[T  =  -  f  p  r2  dv  =  -  f  p  r2 


dv  -  I 


f  .? 


dv 


(Cl  3) 


In  these  integrals,  VR  represents  an  integral  over  any  rigid  parts  of  the 
body  and  the  summation  is  over  all  the  elements.  The  position  vector  of  a 
general  point  in  the  cross  section  of  the  ith  element  can  be  written  as 
follows: 


r  =  °r  +  y°a_2  +  z°aj 


(C14) 


where  °jr  =  °rQ  +  °u_0  (Cl  5) 

As  before,  represents  the  undeformed  position  vector  of  the  cen¬ 
terline  of  the  beam  element  and  °jjq  is  the  initial  displacement  vector  of 
the  centerline.  The  A£  vector  is  neglected  in  these  integrals  since  it  is 
small  and  the  terms  are  of  second  order.  Using  eqs.  (3.12)  and  (3.14),  °r 
is  written  as  follows: 

3 

°r  =  E  Ni(5)(°rQ  +  °uQ)i  (C16) 

“  j=l  J  -0  -0  j 


where 


°uT 


=  [°d  °a 

-oj  H°6j-5  4°6j-4 


>q  1 


(C17) 


Note  that  are  the  x^ ,  X£,  and  x^  displacements  of  the  jth  node. 
Similar  expressions  exist  for  °a^  and  °a^  as  follows: 

3 

°a'  E  MO  °a'  (C18) 

-4  j=1  J  O 

°l3s  X  "J10  '***  (C19) 


The  "al,  and  °al,  terms  can  be  evaluated  using  eqs.  (A7)  and  (A8). 


Thus,  the  finite  element  approximation  to  r.  is  found  by  substituting  eqs. 


(C16)-(C19)  into  eq.  (C14): 


L=  1  Nj(5)  (°r0  +  %  +  y°a^  +  2°a^). 

J  ^ 


(C20) 


If  eq.  (C14)  is  substituted  into  eq.  ( C 1 3)  and  the  area  integrals  carried 


out,  then 


It  =1r  *  °r  d*  +  p 


J  (“ij)2  «  *  »  lyy i]  <^3>2  d'l 


_i  (C21) 


The  matrix  is  the  inertia  matrix  for  the  rigid  parts  of  the  structure 
which  is  found  in  the  usual  way.  The  integrals  along  the  length  of  the 
beam  in  eq.  (C21)  are  evaluated  numerically.  They  can  be  done  exactly 
using  three  point  Gaussian  integration: 


— ?  ^  —2 

°L  d  i  =-2  i  °r.2Uk)  Hk 


(C22) 


l.  3 

1  ^  0-1  2 


°-r  d*  =  T  kf1  °±'  (ck>  Hk 


r  =  2,3 


(C23) 


The  °£(€k)  and  °a_,2(sk)  matrices  can  be  evaluated  by  using  eqs.  (C16)  - 
(C19)  with  £  =  £k. 

The  Element  H  Matrix  (18x3) 


H T  -  /  p  r  N*  dV 


(C24) 
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Substitute  eqs.  (3.18)  and  (C14)  into  (C24)  and  integrate  to  obtain  the 
following  for  the  ith  element: 


H‘  = 


P  M  °il2N$d*+  P  I2ZJ  °£2  dt+  P  Iyy1  /  °a^N*d* 


Li 


Li 


(C25) 


Once  again  eq.  (C25)  is  evaluated  numerically  as  follows: 

j,  rvl«*)N;uk)  ♦  *  Iyy1“i3(«k)il?<5k»Hk 

(C26) 


The  H  matrix  is  assembled  in  the  usual  way, 


The  Mcent  Vector  (3x1 ) 


-cent 


p  r  fi2  r  dv  = 


p  r  a2  r  dv  +  i  I  p  r  ^  r  dv  (C27) 


N 


—2 


i=l 


vi 


Let  I!centR  =  /  p  L  F  X  dv 


'R 


(C28) 


It  can  be  shown  [34]  that  the  eq.  (C28)  is  equivalent  to  the  following: 


(McentRh  =  Ir12  nln3  "  1 R 1 3  nlfi2  “  Ir23  ^2’^  “  ^Ir22_Ir33^  n2fi3 

(C29) 

(McentR^2  =  Ir23  nlQ2  '  Ir12  Q2Q 3  “  IrJ3  (fi3‘^  '  ^Ir33-Ir11^  %nl 

(C30) 
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(McentR^3  =  !R]  3  fi2n3  "  ^23  nl  n3  -  XRi  2  ^n?“^  "  ^IRn‘IR22^  nln2 

(C31 ) 

Note  that  ID  are  the  inertias  of  the  rigid  parts  of  the  structure.  The 
KiJ 

rest  of  Mcent  is  found  as  in  the  previous  sections,  i.e.,  substituting  eq. 
(C 1 4)  into  eq.  (C27)  and  integrating.  Thus, 


N  pi.  3 

Scent  '  ScentR  +  ,f,  ~T~  [Ai  °-S>(5k 1  ^  ”-(Ek' 


*  ‘zz,  “^(EkK^(ek)  *  lyy.  "aJ(Ek)£2"a^({k)]HRl  (C32) 


The  Mcop  Vector  (3x1 ) 


^cor 


N 

£  2 
i  =  l 


p  r  o  u  dv 


(C33) 


Within  an  element,  u  can  be  found  using  eq.  (2.14) 


where 


-y 


“z 


If  a  =  (N£  +  y  N*  +  z  N£)  q  =  +  yUy  +  zuz 

(C34) 

• 

• 

-  N*  q 
—o 

(C35) 

• 

=  N*  q 

(C36) 

-y  - 

• 

• 

=  N*  q 

(C37 ) 
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Substituting  eqs.  (C35)  -  (C37)  as  well  as  eq.  (C14)  into  eq.  (C33)  an 


integrating  results  in  the  following  expression  for  M 


cor 


N  3 

-cor  '  ,f,  »li  (J,  Cfli  W  i£!o<«k> 


+  iZZi  “i2(5kHuk)  *  ‘yy,  °«3<ek)i2i2(C:k)3Hkf 


(C38) 


APPENDIX  D 


DETAILS  OF  THE  MATRICES  AND  VECTORS  RESULTING 

THE  PRINCIPLE  OF  VIRTUAL 


FROM  THE  INERTIA  TERM  IN 
WORK 


The  Element  Mass  Matrix  M  (18x18) 


M  =  /  p  N*  N*  dv 
'V, 


(Dl) 


Substitute  eq.  (3.18)  into  eq.  (Dl),  evaluate  the  matrix  products  and 
integrate  over  the  area  to  obtain 


M  =  Mq  +  +  m  z 


(02) 


where 


Mo  =  p  A  /  N*  N£  d£ 


(03) 


4 


-  P  Izz  /  NJT  NJ  dt 


(04) 


Mz=  P  Iyy  I  I*T  N*  dt 
L 


(D5) 


These  integrals  can  be  evaluated  analytically  using  the  expressions  for  N*, 
N*  and  N_*  in  eqs.  (3.19)  -  (3.21).  The  results  are  as  follows: 
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SYMMETRIC 


0  0  0 
4d}d1  0  2 

0 

16 

SYMMETRIC 


0  0  0 

«iii  2  2 

0 

16 

SYMMETRIC 


E 


3x3  Identity  Matrix 
3  x  3  Null  Matri x 


The  Element  Gyroscopic  Matrix  C 


£  ■  zj  p  N.*T  n  £’ 


N*  dv 


(D9) 


Using  es.  (3.18)  in  Eq.  (09)  and  evaluating  the  integral  as  above  results  in 


C  =  C  +  c  +  C 
—  —o  — y  —  z 


r  -  pAt 

-o  TF~ 


C..  = 


pi  A 


zz 


-y 


IF 


pi  l 


4n 


0  2£ 

£  0 
TBIT 


£ 

£ 

£ 

0 


SYMMETRIC 


■£ 

£ 

If 

£ 

415 


0 


0 


■*£|fiDj  u 


c£|"  QD^ 


0 


0 


SYMMETRIC 


1 6D^ nD2 


£ 

£ 

£ 

0 

0 


0 


4E-|  nE  ^  £ 


2  E 1 9E2 


0 


0 


SYMMETRIC 


1 6E^nEo 


0 

0 

£ 

0 


£ 

£ 

£ 

£ 

£ 

0 


0 


-D[fiD3 


0 


2oI| 03 


4D3OD3 


■l|^3 


0 


2eJoE3 


(010) 


(Dll) 


(012) 


(013) 


The  element  Force  Vector  Due  to  Centrifugal  Force  _FC 


p  N*  £2  (rQ  +  uQ)  dv 


(D14) 


But  by  using  eq.  (C14)  and  (C15), 


I0  +  Hq  =  °£0  +  °M0  +  y°i2  +  2°— 3 


(D15) 


As  in  Appendix  C,  use  numerical  integration  to  evaluate  eq.  ( Dl 4) : 


£c  ■  ^  s  t*  !5T(«k)  i  °K«k)  +  ’zz  i?T('k)  £2  ”^(«k) 


+  ‘yyJSTUk)?  °i3(Ek)]  Hk 


(D16) 
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